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Executive  Summary 


SQUID  NDE  AND  POD  USING  A 
BEM  MEASUREMENT  MODEL 

by  Anthony  P.  Ewing 


Overview 

As  the  commercial  and  military  aircraft  fleets  age,  additional  resources  are  required 
to  ensure  their  airworthiness.  As  the  aircraft  become  older,  the  more  likely  they  are  to 
develop  structural  damage  that  may  lead  to  unscheduled  repairs  or,  in  the  worst  case, 
accidents.  Fatigue  and  corrosion  are  the  two  main  causes  of  structural  damage  in  aging 
aircraft  and  this  research  examines  the  use  of  a  Superconducting  QUantum  Interference 
Device  (SQUID)  magnetometer  as  a  tool  for  Nondestructive  Evaluation  (NDE)  to  detect 
and  characterize  these  aging  aircraft  problems.  The  primary  advantage  of  using  SQUIDs 
in  NDE  over  other  techniques  is  the  ability  to  detect  second  layer  cracks  and  corrosion 
commonly  found  in  aircraft  structures. 

In  general,  verification  of  a  NDE  method  means  demonstrating,  through 
experiment  and/or  calculations,  the  ability  to  distinguish  signal  from  noise  for  the  flaw 
types  and  sizes  and  instrument/flaw  configurations  expected  in  the  actual  inspection.  A 
common  approach  to  quantify  and  validate  the  capabilities  of  an  inspection  technique  is 
to  conduct  a  probability  of  detection  (POD)  analysis.  There  are  basically  two  ways  to 
conduct  this  t5q)e  of  analysis.  The  first  is,  experimentally,  which  requires  a  large  number 
of  samples  with  ranging  flaw  characteristics  being  examined  by  several  inspectors.  The 
second  is,  analytically,  which  requires  a  model  simulating  the  inspection  process  for  a 
range  of  samples  and  testing  conditions. 

A  POD  analysis  has  been  done  using  the  analytical  approach,  combined  with 
experimental  information,  to  evaluate  SQUID  NDE  reliability  consistent  with  damage 
tolerant  design  philosophies  used  in  aerospace  to  make  life  predictions.  A  minimum 
detectable  crack  length  at  90%  probability  of  detection  and  95%  confidence  was  used  as 
the  reliability  criteria. 

Fatigue  Cracks  -  a  Problem  Requiring  NDE 

Fatigue  cracking  in  aircraft  is  primarily  due  to  cyclic  stress  loading.  For  example, 
internal  pressurization  during  flight  creates  stresses  on  the  fuselage  longitudinal  skin 
splices.  The  two  basic  splice  designs  used  are  shown  in  Figure  1 . 
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Fatigue  cracks  are  most  likely  to  develop  along  the  direction  of  the  fastener  row  since  the 
primary  stress  direction  is  transverse  to  the  row  (see  Figure  2).  Usually,  the  cracks  start  at 
the  fastener  hole  and  propagate  radially  and,  if  there  is  multiple  site  damage  (MSD),  the 
structural  strength  is  greatly  reduced  even  for  small  crack  lengths. 
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Figure  2  Fatigue  crack  development  along  direction  of  fastener  row 

Aircraft  are  designed  to  meet  damage  tolerant  requirements  and  nondestructive 
evaluation/inspection  (NDE/I)  techniques  are  needed  to  efficiently  and  reliably  detect  and 
characterize  damage,  if  it  occurs,  so  that  repairs  can  be  made.  Damage  tolerant  design 
basically  involves  fatigue  crack  growth  predictions  starting  from  some  assumed  crack 
length,  Oq  ,  due  to  possible  manufacturing  defects.  For  aircraft  structures,  damage 
tolerance  specifications  require  a  90%  detection  probability  and  95%  confidence  level  for 
detecting  a  specified  crack  length,  (threshold  crack  length),  at  a  particular  location. 

The  dependence  of  a,,,  on  the  inspection  process  and  the  imcertainties  associated  with  the 
detection  of  a  crack  lead  to  the  statistical  problem  that  can  be  addressed  by  POD  analyses. 
Earlier  NDE  work  established  POD  as  a  way  in  which  NDE  process  performance  could 
be  quantified  and  incorporated  into  specifications,  standards  and  design  documents.  By 
quantifying  the  procedure  to  measure  the  performance  capability  of  a  NDE  system, 
objective  comparison  could  be  done  for  different  NDE  systems  and  NDE  performance 
requirements  could  be  defined  for  development  of  new  techniques. 

Probability  of  Detection  (POD) 

Probability  of  detection  can  be  defined  as  the  probability  that  a  specific  crack  length 
can  be  detected  with  a  particular  inspection  system  under  known  conditions.  POD 
combines  characteristics  of  the  measurement  system,  including  noise,  with  statistical 
information  pertinent  to  the  cracks  being  examined.  The  schematic  POD  curve  in  Figure  3 
shows  the  probability  of  detection  as  a  function  of  crack  length.  The  probability  of 
detecting  a  small  crack  is  low  whereas,  the  probability  of  detecting  larger  cracks  approaches 
unity. 
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Crack  Length 

Figure  3  Probability  of  detection  as  a  function  of  crack  length 


The  POD  function  serves  as  a  basis  to  evaluate  the  capability  of  an  NDE  system. 
For  these  analyses,  POD  is  a  function  of  crack  length  to  be  consistent  with  life  prediction 
calculations  that  are  based  on  crack  length.  The  POD  then  is  a  measure  of  how  well  cracks 
of  all  sizes  can  be  found.  This  capability  cannot  be  determined  exactly  but,  through  the  use 
of  confidence  bounds  that  account  for  the  random  errors  of  the  NDE  process,  can  be 
estimated.  For  example,  the  Air  Force  uses  a  95%  lower  confidence  bound  that  provides  a 
risk  factor  that  the  true  probability  of  detection  is  better  than  this  bound  95%  of  the  time. 


Use  of  SQUID’S  in  NDE 

A  superconducting  quantum  interference  device  (SQUID)  magnetometer  is  the 
world’s  most  sensitive  detector  of  magnetic  fields.  The  potential  advantage  that  a  SQUID 
has  over  other  techniques  is  its  ability  to  detect  second  layer  flaws  by  using  low 
frequency  excitation. 

A  SQUID  NDE  technique  is  not  only  defined  by  the  SQUID  instrument 
characteristics  but  also  by  nature  of  the  current  distributions  used  to  probe  the  flaws. 
When  the  current  passing  through  the  sample  is  affected  by  inhomogenieties,  voids, 
cracks,  and  edges,  the  magnetic  field  is  perturbed  and  can  be  sensed  by  the  pickup  coil  of 
the  SQUID.  This  work  focuses  on  the  dc-current  injection  technique.  Injecting  or  inducing 
a  uniform  current  distribution  in  the  specimen  causes  the  current  to  be  parallel  to  the 
specimen  surface  under  the  pickup  coil.  For  example,  a  large,  uniform  plate  would  have 
a  magnetic  field  from  this  that  encircles  the  specimen.  This  magnetic  field  would 
primarily  be  parallel  to  the  specimen  surface  for  scans  centrally  located  and  for  small 
liftoff  distances.  The  pickup  coils  measure  only  the  perpendicular  component  (or 
gradient  of  this)  of  the  magnetic  field.  A  flaw  in  the  specimen  will  perturb  the  parallel 
field  and  produce  a  perpendicular  component  which  can  then  be  detected  (pickup  coils 
measure  only  BJ.  Figure  4  illustrates  the  SQUID  pickup  coil  being  scanned  over  the 
sample  containing  a  flaw  and  the  typical  magnetic  map  produced  revealing  a  signature 
that  commonly  has  a  dipolar  shape.  The  trough  and  ridge  at  the  ends  of  the  map  are  due 
to  the  edge  effects  of  the  plate. 
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Figure  4  Dc-current  injection:  (a)  scan  over  sample  (b)  resulting  magnetic  map 

Although  this  technique  is  not  likely  to  be  used  in  a  field  instrument,  the  current 
distributions  produced  are  similar  to  the  planar  ac  eddy  currents  produced  by  sheet  inducers, 
which  are  proposed  for  a  field  instrument.  Also,  since  existing  boundary  element  methods 
can  be  used  to  model  dc-current  injection,  the  analytical  approach  to  a  POD  analysis  can  be 
taken. 

Measurement  Model 

A  measurement  model  using  boundary  element  methods  has  been  constructed 
simulating  a  SQUID  magnetometer  being  scanned  over  a  dc-current  injected,  finite  plate 
containing  an  ideal  {i.e.,  straight,  infinitely  thin,  and  perfectly  insulated)  crack.  The 
model  was  used  to  examine  the  effect  of  system  parameter  variability  on  flaw 
detectability.  The  results  of  this  sensitivity  analysis  was  then  combined  with  empirical 
noise  distributions  to  determine  probability  of  detection.  A  goal  in  this  research  was  to 
develop  a  method  where  simulation  represents  the  experimental  approach,  is  cheaper  and 
faster,  and  identifies  sources  of  unreliability  in  SQUID  NDE. 

The  objective  of  the  SQUID  measurement  model,  which  solves  the  forward 
problem  of  calculating  the  magnetic  field  from  a  known  current  distribution,  is  to 
simulate  what  the  instrument  may  see  in  a  test  environment  on  an  unknown  sample.  A 
SQUID  gradiometer  measures  the  vertical  component  of  magnetic  field,  which  is  the 
result  of  currents  flowing  in  the  sample.  A  flaw  will  perturb  these  currents  and,  if  the 
perturbation  is  large  enough  {e.g.,  signal/noise  >1),  show  up  as  an  anomaly  in  the 
magnetic  field  map.  The  specific  aging  aircraft  problem  being  addressed  is  the  modeling 
of  a  fatigue  crack  emanating  radially  from  a  fastener  hole  (schematically  shown  in  Fig.  5) 
as  discussed  previously  regarding  fuselage  skin  splices. 


Fastener  hole 


Figure  5  Schematic  of  a  fastener  hole  with  crack 


For  the  interested  reader,  the  details  of  the  mathematical  development  of  this  can 
be  found  in  Chapter  II  of  the  associated  dissertation  from  which  this  summary  is  written. 

POD  Curves  from  Simulated  Experimental  Data 

Figure  6  shows  the  centerline  profiles  of  the  magnetic  maps  resulting  from 
simulated  scans  over  a  hole  with  crack  for  four  different  crack  lengths  while  keeping  the 
hole  diameter  constant.  As  can  be  seen,  for  increasing  crack  length,  the  signal 


amplitude  increases.  Also,  the  asymmetry  increases,  pulling  the  crack  side  of  the  dipolar 
signal  to  the  right  while  the  hole  side  remains  relatively  fixed.  The  important  characteristic 
for  POD  is  the  increasing  peak-to-peak  value  as  a  fimction  of  crack  length,  which  is  the 
“signal”  used  in  the  POD  analysis. 

Sensitivity  analyses  were  done  on  the  system  parameters  of  scan  resolution,  plate 
thickness,  pick-up  coil  liftoff,  and  pick-up  coil  tilt  angle  to  determine  how  they  affect  the 
overall  signal  distribution.  Monte  Carlo  simulations  used  this  signal  distribution  along 
'with  experimental  noise  distributions  to  generate  realistic  POD  curves. 

Noise  distributions  must  be  characterized  in  order  to  determine  POD  and  POFA. 
Several  noise  distributions  have  been  extracted  from  existing  experimental  data  for  different 
SQUID  systems  and  measurement  techniques.  The  measurement  model  simulates  dc- 
current  injection  and  experiments  using  direct  dc-current  injection  have  large  noise 
distributions  since  they  use  the  entire  bandwidth  of  the  SQUID  (dc  to  -lOkHz)  and 
therefore,  include  noise  over  these  frequencies  as  well.  Noise  conditions  associated  with 
other  experimental  techniques  are  more  representative  of  what  SQUID’s  will  be  operated  in. 
Techniques  based  upon  eddy  current  inducers  use  lock-in  amplifiers  at  a  particular 
frequency  and  greatly  reduce  the  noise.  The  dc-measurement  data  shows  approximately 
400  pT  peak-to-peak  noise  while  the  lock-in  measurement  data  shows  about  2  pT,  a  noise 
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reduction  factor  of  200!  POD  values  determined  using  noise  from  lock-in  measurements 
represent  present  SQUID  capability. 

Monte  Carlo  simulation,  utilizing  the  results  from  the  BEM  measurement  model, 
was  used  to  sample  from  the  imcertainty  distributions  of  scan  resolution,  plate  thickness, 
pick-up  coil  liftoff  and  pick-up  coil  tilt  angle  to  generate  an  overall  signal  distribution.  The 
resulting  signal  distribution  was  then  compared  with  the  noise  distributions  associated  with 
ac  and  dc  measurements  to  determine  POD.  A  95%  lower  confidence  bound  was  then 
generated  by  iterative  generation  of  multiple  signal  distributions  at  each  crack  length. 

Figure  7  is  the  dc-measurement  POD  curve  resulting  from  generating  similar 
distributions  at  each  crack  length.  The  abruptness  of  the  POD  curve  near  the  origin  is  due 
to  a  combination  of  the  relative  sharpness  of  the  signal  distribution  with  respect  to  the  noise 
distribution  and  the  one-sidedness  of  the  noise  distribution.  The  95%  lower  confidence 
limit  applies  more  to  experimentally  based  POD  curves  that  are  usually  generated  from  a 
relatively  small  number  of  data  points.  For  simulated  data,  the  lower  confidence  bound 
can  be  made  to  basically  lie  on  the  POD  curve  if  enough  runs  are  made  and  is  not  as 
useful  as  a  concept  as  for  the  experimentally  derived  POD  curve. 


Crack  Length  (mm) 

Figure  7  Ac/dc-measurement  POD  curves  showing  90/95  crack  lengths 
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The  minimum  detectable  crack  lengths  corresponding  to  90%  probability  of 
detection  at  95%  confidence  are  1.4  mm  for  the  dc-measurement  and  0.0134  mm  for  the  ac- 
measurement.  The  very  small  minimum  detectable  crack  length  determined  for  the  ac- 
measurement  is  due  to  the  large  noise  reduction  through  use  of  the  lock-in  amplifier  in  this 
type  of  measurement. 

In  the  Monte  Carlo  simulations,  the  signal  distribution  was  compared  directly  with 
the  noise  distribution  to  determine  POD.  For  direct  sampling  from  distributions,  the  signal- 
to-noise  ratio  set  the  threshold  ( a,^  ),  which  set  the  detection  criteria.  To  examine  the 
tradeoff  between  the  probability  of  detection  (POD)  and  the  probability  of  false  alarm 
(POFA)  when  setting  ,  the  signal  and  noise  distribution  data  generated  by  the  Monte 
Carlo  simulations  can  be  presented  in  a  different  format.  The  SNR  requirements  determine 
where  a,^  is  to  be  placed,  thus  affecting  POD  (larger  SNR  requirements  correspond  to  larger 
minimum  detectable  crack  lengths). 

Figure  8  displays  the  dc-measurement  POFA  and  POD  curves  for  three  crack 
lengths,  including  the  minimum  detectable  crack  length,  as  a  fimction  of  threshold.  By 
setting  a  threshold,  a  value  for  POFA  and  POD  for  all  crack  lengths  is  determined.  As 
can  be  seen,  for  a  crack  length  of  0.5  mm,  the  90%  probability  of  detection  point 
corresponds  to  a  28%  POFA  (i.e.,  there  is  a  28%  chance  that  noise  will  be  mistaken  for  a 
signal).  There  is  only  a  single  POFA  curve  since  probability  of  false  alarm  is  determined 
from  the  noise  distribution,  which  for  these  simulations  is  constant.  Figure  9  shows  the 
ac-measurement  POFA  and  POD  as  a  function  of  threshold  for  the  previously  determined 
minimum  detectable  crack  length  of  0.0134  mm.  In  both  the  dc  and  ac  measurement 
cases,  the  POFA  corresponding  to  the  minimum  detectable  crack  length  is  approximately 
zero. 


Threshold  (nT) 

Figure  8  Dc-measurement  POFA  and  POD  as  a  function  of  threshold 
for  three  values  of  crack  length 
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Figure  9  Ac-measurement  POFA  and  POD  as  a  function  of  threshold 
for  the  minimum  detectable  crack  length 


Factors  That  May  Decrease  POD  Capability 

The  relatively  small  values  of  minimum  detectable  crack  length  determined  from 
this  analysis  prompts  discussion  on  whether  this  reflects  the  true  capability  of  SQUIDs  in 
the  detection  of  real  cracks  near  fasteners  in  aircraft  lap  joints.  Experimental  measurements 
using  SQUIDs  on  fatigue  cracks  in  lap  joints  have  demonstrated  the  difficulty  in  finding 
cracks,  adjacent  to  fasteners,  on  the  order  several  millimeters  in  length.  For  the  POD 
analysis  done  here,  only  four  system  parameters  (scan  resolution,  plate  thickness,  pickup 
coil  liftoff  and  tilt  angle)  and  their  associated  uncertainties  were  considered.  Since  it  is 
these  parameter  uncertainties  that  determine  the  overall  signal  distribution,  and  therefore 
POD,  it  is  apparent  that  significant  parameters  have  been  left  out  of  this  analysis.  More 
than  likely,  these  stem  from  sample-related  noise  associated  with  real  measurements,  which 
have  not  been  represented  in  either  the  measurement  model  itself  or  through  the  parameter 
uncertainty  distributions.  Sample-related  noise  results  from  the  geometry  of  the  sample  or, 
more  specifically,  how  that  geometry  affects  the  current  flow  used  to  probe  the  sample. 
Figure  10  shows  how  some  of  these  current  paths  can  go  under  and/or  through  the  crack  and 
the  fastener,  further  complicating  the  detection  and  characterization  of  the  crack.  Some  of 
these  issues  are  discussed  in  the  following  sections. 
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Figure  10  Additional  current  paths  complicating  crack  detection 


Crack  shape  -  Probability  of  detection  is  traditionally  stated  as  a  function  of  crack 
length  to  be  consistent  with  damage  tolerant  requirements,  which  are  stated  in  terms  of 
crack  length.  However,  crack  shape  is  an  important  issue  to  discuss  briefly  since  cracks  of 
the  same  “length”  but  different  shape  could  result  in  different  signals.  Currently,  the  BEM 
measurement  model  simulates  a  straight-through  crack  (the  same  length  through  the  entire 
thickness  of  the  plate)  but  it  is  possible  to  use  the  model  to  simulate  “3-D”  cracks  to 
examine  the  effect  of  crack  “shape”  on  magnetic  field.  By  creating  slices  of  varying  shapes 
and  depths  (see  Fig.  1 1),  a  pseudo-3D  crack  can  be  modeled. 


Figure  1 1  Pseudo  3-D  modeling  of  a  crack  using  slices 

Since  crack  shape  affects  both  the  amplitude  and  the  shape  of  the  magnetic  field,  then  this 
pseudo  3-D  approach  would  be  a  better  approximation  to  a  real  crack  that  has  shape  that 
varies  with  depth. 

Ideal  cracks  versus  real  cracks  and  slots  -  The  BEM  measurement  model  assumes 
an  ideal,  perfectly  insulating,  and  infinitely  thin  crack.  This  assumption  will  lead  to 
differences  between  the  model  and  measurements  made  on  samples  containing  either 
wide  cracks  (slots)  or  real  fatigue  cracks.  Most  fabricated  test  samples  used  in  our 
laboratory  use  slots,  made  by  a  saw  or  electrodischarge  machining  (EDM),  to  simulate 
cracks.  A  combination  of  a  drilled  hole  with  an  EDM  slot  is  an  approximation  to  a  crack 
emanating  from  an  aircraft  fastener  hole.  This  is  useful  in  preliminary  analyses,  especially 
for  development  of  NDE  techniques.  Fabrication  of  test  samples  made  this  way  is 
simple  and  controllable,  making  it  easy  to  build  a  test  set  representing  the  range  of 
conditions  that  are  of  interest.  But,  measurements  with  NDE  instruments  have  shown 
that  the  instrument  response  from  a  slot  is  not  necessarily  the  same  as  that  from  a  fatigue 
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crack  of  the  same  size  and  geometry.  Figure  12  schematically  shows  some  of  the 
possible  variations,  each  possibly  resulting  in  a  different  instrument  response. 


Figure  12  Slot  profiles:  (a)  wide  slot  (b)  narrow  slot  (c)  ideal  crack  (d)  real  crack 


The  BEM  measurement  model  simulates  a  closed  crack  that  is  electrically  insulated 
along  the  entire  length.  We  have  shown  that  it  provides  a  reasonable  approximation  to  a 
thin  slot  for  liftoffs  appropriate  for  SQUIDs  commonly  used  in  the  laboratory  (>  3  mm). 
However,  for  real  fatigue  cracks,  it  is  possible  that  crack  closure  may  cause  electrical 
conductivity  across  parts  of  the  crack.  Probability  of  detection  will  be  strongly  dependent 
on  the  effects  of  crack  closure  since,  if  current  is  flowing  through  the  crack  instead  of 
around  it,  the  signal  (which  is  proportional  to  the  perturbation  of  the  current)  will  be  greatly 
reduced.  It  is  not  yet  understood  how  much  crack  closure  effects  current  flow  since  it  is 
likely  that  the  oxide  layers  that  form  on  the  crack  surfaces  will  act  as  an  electrical  insulator 
and  so,  only  those  contact  points  where  the  oxide  layer  is  not  present  will  provide 
conductive  paths. 

Additional  current  paths  -  The  POD  analysis  has  been  based  on  simulated 
measurements  of  a  crack  emanating  from  a  fastener  hole  and  has  completely  ignored  the 
effects  of  the  fastener  itself  This  is  most  likely  one  of  the  larger  sources  of  discrepancy 
present  between  the  model  and  real  measurement.  Since  most  fatigue  cracks  start  imder  the 
fastener  head,  early  detection  is  difficult  and  usually  the  crack  has  to  propagate  beyond  the 
fastener  head  before  detection  occurs.  Another  critical  issue  is  that  of  contact  between  the 
fastener  and  the  hole  sides,  which  can  create  numerous  current  flow  paths  across  and 
through  the  fastener.  These  currents  may  be  too  difficult  to  model  directly  and  their 
associated  uncertainty  might  have  to  be  represented  through  a  sample-related  noise 
distribution.  This  will  probably  have  to  be  determined  experimentally  by  examining 
samples  with  holes  alone  and  with  fasteners  inserted  to  quantify  the  effect  on  the  magnetic 
signal. 

Geometry  factors  -  Current  can  also  be  greatly  affected  by  the  geometry  of  the 
surrounding  structure,  particularly  edges  such  as  lap  joint  seams.  Figure  13  shows  that 
edges  can  produce  a  large  signal  amplitude  that  can  make  flaw  detection  difficult,  especially 
if  Ae  flaw  is  located  near  an  edge.  Image  processing  techniques  can  partially  remove  the 
background  slope  due  to  these  edges  but  it  is  still  difficult  to  extract  the  signal  from  a  flaw 
that  is  near  an  edge.  Edges  and  other  geometrical  factors  {e.g.,  structural  support  members) 
affect  the  ability  of  the  SQUID  to  detect  a  crack  and  these  have  not  been  taken  into  account 
in  this  POD  analysis. 
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Figure  13  Magnetic  map  centerline  profile  pointing  out  edge  and  flaw  signals 

Factors  That  May  Increase  POD  Capability 

Signal  definition  -  The  peak-to-peak  value  of  the  magnetic  dipolar  map  is  the 
present  definition  of  “signal”  used  in  these  analyses.  Peak-to-peak  amplitude  is  a  very 
limited  use  of  the  information  available  in  the  mapping  since  it  reduces  the  entire  2-D  image 
to  a  single  value.  It  is  the  opinion  of  the  author  that  image  features  {e.g. ,  shape  and 
asymmetry)  of  the  2-D  magnetic  field  map  are  just  as  important  as  the  peak-to-peak  signal 
amplitude  in  the  detection  and  characterization  of  cracks,  especially  near  fastener  holes.  For 
example,  if  we  look  at  the  magnetic  map  centerline  profiles  for  various  ideal  geometries  that 
have  the  same  “cross  sectional  length”  perpendicular  to  the  current  injection  direction,  we 
can  see  that  using  only  peak-to-peak  information  for  current  perpendicular  to  the  crack  can 
lead  to  detection  problems.  The  left-hand  side  of  Fig.  14  shows  the  profiles  for  the  current 
transverse  to  the  crack  at  three  different  liftoffs  for  a  9  mm  diameter  hole,  a  5  mm  hole  with 
a  4  mm  crack,  and  a  9  mm  crack.  The  signal  values  for  each  of  the  traces  are  plotted  on  the 
same  scale  and  represent  a  SQUID  system  using  3  mm  pickup/balance  coils  with  a  3  cm 
baseline.  From  this  information,  all  the  profiles  look  the  same  in  that  they  are  dipolar  with  a 
mostly  symmetrical  shape.  Based  on  peak-to-peak  amplitude  in  a  single  direction,  we  can 
not  distinguish  between  a  crack  alone,  hole  alone,  or  a  hole-with-crack  combination  imless  a 
baseline  signal  amplitude  can  be  established  as  was  done  earlier  in  this  POD  analysis  using 
the  signal  associated  with  a  5  mm  hole  alone.  However,  this  approach  would  not  work  here 
since  the  5  mm  hole  is  not  a  constant  feature  in  all  three  geometries. 

The  rotating  current  schemes  take  advantage  of  the  differing  two-dimensional 
structme  of  these  three  geometries.  By  rotating  the  current  direction,  the  signals  associated 
with  a  crack  will  go  through  a  maximum  (when  the  current  is  perpendicular  to  the  crack) 
and  a  minimum  (when  the  current  is  parallel  to  the  crack).  A  hole  alone  will  not  show  this 
cyclic  behavior  (unless  it  is  out-of-round).  On  the  right-hand  side  of  Fig.  14  are  the 
corresponding  profiles  when  the  current  is  parallel  to  the  crack  showing  the  9  mm  crack 
having  a  flatline  (zero  signal)  and  the  5  mm  hole  with  a  4  mm  crack  having  a  signal 
corresponding  to  a  5  mm  hole  alone  and  the  9  mm  hole  is  unchanged.  In  this  way,  the  three 
geometries  could  be  distinguished  by  a  measure  of  their  peak-to-peak  difference  at  the  two 
current  injection  directions,  with  the  crack  alone  being  the  largest,  the  hole-with-crack 
combination  the  next  largest,  and  the  hole  alone  having  no  difference.  Real  fatigue  cracks 
will  have  some  small  signal  for  current  injected  parallel  to  the  crack  (versus  zero  signal  for 
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Signal  (arbitrary  units) 


the  ideal  cracks  of  Fig.  14)  but  will  still  show  the  same  maximum/minimum  behavior  as  a 
function  of  current  injection  direction  (the  signal  difference  will  just  be  smaller  than  the 
ideal  crack  case). 


Figure  14  Profiles  for  three  9-mm  “cross  section”  geometries  at  varying  liftoffs 

(Note:  all  profile  traces  are  on  the  same  scale) 

SQUID  NDE  systems  will  have  to  use  asymmetry  measures  (which  includes  the 
rotating  current  direction  technique)  for  crack  detection,  especially  for  cracks  originating 
from  fastener  holes.  If  the  asymmetry  technique  is  used  only  to  find  the  maximum  and 
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minimum  of  the  signal,  then  POD  may  still  be  determined  primarily  from  peak-to-peak 
amplitudes.  However,  if  additional  measures  of  asymmetry  {e.g.,  eccentricity  or  other 
shape  factors)  are  to  be  used,  then  POD  analyses  must  incorporate  these  measures  into  the 
standard  peak-to-peak  analysis  to  reflect  the  true  capability  of  the  system.  This  is  an 
important  topic  to  be  addressed  for  fiiture  development  of  SQUID  POD  methodologies. 

Conclusions 

A  measurement  model  using  boundary  element  methods  has  been  constructed 
simulating  a  SQUID  magnetometer  being  scanned  over  a  dc-current  injected,  finite  plate 
containing  an  ideal  {i.e.,  straight,  infinitely  thin,  and  perfectly  insulated)  crack.  The 
model  was  used  to  examine  the  effect  of  system  parameter  (scan  resolution,  plate 
thickness,  pick-up  coil  liftoff,  and  pick-up  coil  tilt  angle)  variability  on  flaw  detection. 
The  results  of  this  sensitivity  analysis  were  then  combined  with  empirical  noise 
distributions  to  determine  probability  of  detection.  A  goal  in  this  research  was  to  develop 
a  method  where  simulation  represents  the  experimental  approach,  is  cheaper  and  faster, 
and  identifies  sources  of  unreliability  in  SQUID  NDE. 

This  POD  analysis  resulted  in  minimum  detectable  crack  lengths  corresponding  to 
90%  probability  of  detection  at  95%  confidence  of  1 .4  mm  for  the  dc-measurement  and 
0.0134  mm  for  the  ac-measurement.  The  very  small  minimum  detectable  crack  length 
determined  for  the  ac-measurement  is  due  to  the  large  noise  reduction  through  use  of  the 
lock-in  amplifier  in  this  type  of  measurement. .  However,  there  is  a  large  discrepancy 
between  experimental  measurements  on  real  fatigue  cracks  and  the  results  of  this  analysis. 
If  the  model  is  to  be  used  to  calculate  realistic  probability  of  detection  values  for  SQUID 
systems,  more  work  has  to  be  done  to  quantify  additional  sources  of  noise  to  be 
incorporated  into  the  analysis.  Some  of  the  possible  sources  of  this  noise  have  been 
discussed  but  no  work  has  been  done  as  of  yet  to  quantify  their  effects. 

The  BEM  measurement  model  has  already  been  used  in  other  experimental  work 
involving  validation  and  calibration  techniques.  The  model  is  just  now  beginning  to  be 
used  for  POD  work,  which  was  what  it  was  originally  designed  for,  and  continued 
development  will  improve  the  model  as  a  tool  for  SQUID  NDE  development  and 
capability  measurement. 
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Superconducting  QUantum  Interference  Device’s  (SQUID’s)  are  being  used  as  tools 
for  Nondestructive  Evaluation  (NDE)  to  detect  and  characterize  defects  in  aging  aircraft.  To 
evaluate  SQUID  NDE  reliability,  a  probability  of  detection  (POD)  analysis  has  been  done. 

A  boundary  element  method  (BEM)  measurement  model  using  a  Green’s  function 
developed  specifically  for  crack  problems  has  been  constructed  for  use  in  the  POD  analysis. 
The  model  simulates  the  2-D  images  of  the  magnetic  field  obtained  by  scanning  a  SQUID 
magnetometer  over  a  plate  containing  a  crack  and  carrying  an  injected  dc-current.  POD 
curves  were  generated  through  Monte  Carlo  simulation  using  distributions  derived  from 
sensitivity  analyses  and  experimental  noise  measurements.  For  the  conditions  simulated, 
crack  lengths  of  1 .4  mm  (dc-measurement)  and  0.0134  mm  (ac-measurement)  could  be 
found  with  90%  probability  of  detection  and  95%  confidence.  These  small  crack  lengths 
suggest  that  additional  experimental  noise  factors  will  have  be  incorporated  into  the  POD 
analysis  before  realistic  S/^ID  NDE  capability  can  be  accurately  quantified. 
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CHAPTER  I 


INTRODUCTION 


Overview 

As  the  commercial  and  military  aircraft  fleets  age,  additional  resources  are  required 
to  ensure  their  airworthiness.  As  the  aircraft  become  older,  the  more  likely  they  are  to 
develop  structural  damage  that  may  lead  to  unscheduled  repairs  or,  in  the  worst  case, 
accidents.  Fatigue  and  corrosion  are  the  two  main  causes  of  structural  damage  in  aging 
aircraft  and  this  research  examines  the  use  of  a  Superconducting  QUantum  Interference 
Device  (SQUID)  magnetometer  as  a  tool  for  Nondestructive  Evaluation  (NDE)  to  detect 
and  characterize  these  aging  aircraft  problems.  The  primary  advantage  of  using  SQUIDs 
in  NDE  over  other  techniques  is  the  ability  to  detect  second  layer  cracks  and  corrosion 
commonly  foimd  in  aircraft  structures. 

In  general,  verification  of  a  NDE  method  means  demonstrating,  through 
experiment  and/or  calculations,  the  ability  to  distinguish  signal  from  noise  for  the  flaw 
types  and  sizes  and  instrument/flaw  configurations  expected  in  the  actual  inspection.  A 
common  approach  to  quantify  and  validate  the  capabilities  of  an  inspection  technique  is 
to  conduct  a  probability  of  detection  (POD)  analysis.  There  are  basically  two  ways  to 
conduct  this  type  of  analysis.  The  first  is,  experimentally,  which  requires  a  large  number 
of  samples  with  ranging  flaw  characteristics  being  examined  by  several  inspectors.  The 
second  is,  analytically,  which  requires  a  model  simulating  the  inspection  process  for  a 
range  of  samples  and  testing  conditions. 
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A  POD  analysis  is  being  done  using  the  analytical  approach,  combined  with 
experimental  information,  to  evaluate  SQUID  NDE  reliability  consistent  with  damage 
tolerant  design  philosophies  used  in  aerospace  to  make  life  predictions.  A  minimum 
detectable  crack  length  at  90%  probability  of  detection  and  95%  confidence  is  being  used 
as  the  reliability  criteria. 


Fatigue  Cracks  -  a  Problem  Requiring  NDE 

Fatigue  cracking  in  aircraft  is  primarily  due  to  cyclic  stress  loading.  For  example, 
internal  pressurization  during  flight  creates  stresses  on  the  fuselage  longitudinal  skin 
splices.  The  two  basic  splice  designs  used  [1]  are  shown  in  Figure  1.1. 


Figure  1.1  Fuselage  skin  splices:  a)  lap  splice  and  b)  butt  splice 

Fatigue  cracks  are  most  likely  to  develop  along  the  direction  of  the  fastener  row  since  the 
primary  stress  direction  is  transverse  to  the  row  (see  Figure  1 .2).  Usually,  the  cracks  start  at 
the  fastener  hole  and  propagate  radially  and,  if  there  is  multiple  site  damage  (MSD),  the 
structural  strength  is  greatly  reduced  even  for  small  crack  lengths. 
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Figure  1.2  Fatigue  crack  development  along  direction  of  fastener  row 


Aircraft  are  designed  to  meet  damage  tolerant  requirements  [2-5]  and  nondestructive 
evaluation/inspection  (NDE/I)  techniques  are  needed  to  efficiently  and  reliably  detect  and 
characterize  damage,  if  it  occurs,  so  that  repairs  can  be  made.  Damage  tolerant  design 
basically  involves  fatigue  crack  growth  predictions  starting  from  some  assumed  crack 
length,  Oq  ,  due  to  possible  manufacturing  defects.  For  aircraft  structures,  damage 
tolerance  specifications  require  a  90%  detection  probability  and  95%  confidence  level  for 
detecting  a  specified  crack  length,  a,^ ,  at  a  particular  location.  The  inspection  interval,  t, 
is  determined  by 


where  Np  is  the  fracture  mechanics  propagation  life  or  the  life  for  a  crack  length  a  to 
propagate  to  the  permissible  crack  length,  Op,  under  expected  usage  environments  and  S/ 
is  the  safety  factor.  The  permissible  crack  length  Oj ,  which  is  smaller  than  the  failure 
crack  length,  is  associated  with  minimum  residual  strength  under  the  presence  of  cracks, 
defined  by  Federal  Aviation  Regulations  Airworthiness  Requirements  (FAR-25)[6]. 
Commonly,  a  safety  factor  of  two  is  used  and  so  inspections  are  then  required  at  half  the 
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predicted  time  for  the  crack  to  reach  the  permissible  crack  length.  After  the  first 
inspection,  the  crack  length  is  reset  to  the  largest  crack  length  that  can  be  missed  by  the 
inspection  process  being  used,  a,^  (see  Figure  1 .3). 


Figure  1.3  Crack  growth  and  inspection 


This  process  is  repeated  throughout  the  aircraft’s  life  with  detected  cracks  being  repaired. 
The  dependence  of  on  the  inspection  process  and  the  uncertainties  associated  with  the 
detection  of  a  crack  lead  to  the  statistical  problem  that  can  be  addressed  by  POD  analyses. 
Earlier  NDE  work  [7,8]  established  POD  as  a  way  in  which  NDE  process  performance 
could  be  quantified  and  incorporated  into  specifications,  standards  and  design  documents 
[3-5].  By  quantifying  the  procedure  to  measure  the  performance  capability  of  a  NDE 
system,  objective  comparison  could  be  done  for  different  NDE  systems  and  NDE 
performance  requirements  could  be  defined  for  development  of  new  techniques. 
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Probability  of  Detection  (POD) 


Probability  of  detection  can  be  defined  as  the  probability  that  a  specific  crack  length 
can  be  detected  with  a  particular  inspection  system  under  known  conditions.  POD 
combines  characteristics  of  the  measurement  system,  including  noise,  with  statistical 
information  pertinent  to  the  cracks  being  examined.  The  schematic  POD  curve  in  Figure 
1.4  shows  the  probability  of  detection  as  a  function  of  crack  length.  The  probability  of 
detecting  a  small  crack  is  low  whereas,  the  probability  of  detecting  larger  cracks  approaches 
unity. 


Figure  1 .4  Probability  of  detection  as  a  function  of  crack  length 


The  POD  function  serves  as  a  basis  to  evaluate  the  capability  of  an  NDE  system. 
For  these  analyses,  POD  is  a  function  of  crack  length  to  be  consistent  with  life  prediction 
calculations  that  are  based  on  crack  length.  The  POD  then  is  a  measure  of  how  well  cracks 
of  all  sizes  can  be  found.  This  capability  cannot  be  determined  exactly  but,  through  the  use 
of  confidence  bounds  that  accoimt  for  the  random  errors  of  the  NDE  process,  can  be 
estimated.  For  example,  the  Air  Force  uses  a  95%  lower  confidence  bound  that  provides  a 
risk  factor  that  the  true  probability  of  detection  is  better  than  this  bound  95%  of  the  time. 
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Use  of  SQUID’S  in  NDE 


A  superconducting  quantum  interference  device  (SQUID)  magnetometer  is  the 
world’s  most  sensitive  detector  of  magnetic  fields.  The  potential  advantage  that  a  SQUID 
has  over  other  techniques  is  its  ability  to  detect  second  layer  flaws  by  using  low 
frequency  excitation.  Figure  1.5  shows  one  of  the  imaging  SQUID  gradiometer  systems 
at  Vanderbilt  University  (termed  p-SQUID  [9])  which  contains  four  3-mm  diameter,  16 
turn,  closely-spaced  pick-up  coils  for  high  resolution  measurements.  Minimum  liftoff 
(pick-up  coil  to  sample)  distance  is  2-3  mm.  To  reject  noise,  all  pick-up  coils  are  axial 
gradiometers,  each  having  a  counter- wound  12-nim,  1  turn,  balance  coil  with  a  baseline  of 
3  cm.  Two  separate  scanning  stages  are  currently  being  used  to  move  samples  beneath 
SQUID  instruments  to  produce  the  2-D  magnetic  field  maps.  One,  a  non-magnetic  high¬ 
speed  (up  to  10  cm/s)  belt  driven  stage,  is  used  inside  the  magnetic  shield.  The  other,  an 
open  frame  screw-drive  stage,  cannot  be  used  inside  the  shield  and  must  use  a  diving 
board  assembly  to  hold  the  sample  under  the  SQUID  for  scanning. 


Figure  1 .5  The  p-SQUID  system 
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A  SQUID  NDE  technique  is  not  only  defined  by  the  SQUID  instrument 
characteristics  but  also  by  nature  of  the  current  distributions  used  to  probe  the  flaws. 
When  the  current  passing  through  the  sample  is  affected  by  inhomogenieties,  voids, 
cracks,  and  edges,  the  magnetic  field  is  perturbed  and  can  be  sensed  by  the  pickup  coil  of 
the  SQUID.  Some  of  the  commonly  used  techniques  in  SQUID  NDE  include  the 
following: 

a)  Ac/dc  Current  Injection:  Injecting  or  inducing  a  uniform  current  distribution 
in  the  specimen  causes  the  current  to  be  parallel  to  the  specimen  surface  under 
the  pickup  coil.  For  example,  a  large,  uniform  plate  would  have  a  magnetic 
field  from  this  that  encircles  the  specimen.  This  magnetic  field  would 
primarily  be  parallel  to  the  specimen  surface  for  scans  centrally  located  and 
for  small  liftoff  distances.  The  pickup  coils  measure  only  the  perpendicular 
component  (or  gradient  of  this)  of  the  magnetic  field.  A  flaw  in  the  specimen 
will  perturb  the  parallel  field  and  produce  a  perpendicular  component  which 
can  then  be  detected  (pickup  coils  measure  only  BJ.  Figure  1.6  illustrates  the 
SQUID  pickup  coil  being  scanned  over  the  sample  containing  a  flaw  and  the 
typical  magnetic  map  produced  revealing  a  signature  that  commonly  has  a 
dipolar  shape.  The  trough  and  ridge  at  the  ends  of  the  map  are  due  to  the  edge 
effects  of  the  plate. 

b)  Inducing  Plate:  An  eddy-current  technique  using  a  striped  inducer  sheet 
constructed  from  a  300  pm  copper  sheet  carrying  a  low  frequency  ac-current 
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Figure  1.6  Dc-current  injection:  (a)  scan  over  sample  (b)  resulting  magnetic  map 


has  been  used  to  image  flaws  in  test  specimens.  This  also  leads  into  the  present 
work  related  to  the  use  of  phase  information  of  the  eddy  current  image  to 
determine  depth  of  the  flaw  [10]. 

c)  Ac  Nulling:  This  is  basically  a  modified  eddy-current  technique.  An  ac 
magnet  located  in  the  tail  of  the  instrument  is  capable  of  producing  a  low- 
ffequency  ac-signal  (200  Hz  max).  Optimizing  software  controls  feedback 
circuitry  to  cancel  the  ac-signal  that  is  directly  coupled  to  the  SQUID  as  different 
specimens  are  examined.  This  feedback  allows  a  “zeroing”  out  of  the  signal 
before  a  scan  and  enables  the  instrument  to  be  set  on  a  higher  sensitivity. 

General  Approach 

This  work  focuses  on  the  dc-current  injection  technique.  Although  this  technique  is 
not  likely  to  be  used  in  a  field  instrument,  the  current  distributions  produced  are  similar  to 
the  planar  ac  eddy  currents  produced  by  sheet  inducers,  which  are  proposed  for  a  field 
instrument.  Also,  since  existing  boundary  element  methods  can  be  used  to  model  dc- 
current  injection,  the  analytical  approach  to  a  POD  analysis  can  be  taken. 
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A  measurement  model  using  boundary  element  methods  will  be  constructed 
simulating  a  SQUID  magnetometer  being  scanned  over  a  dc-current  injected,  finite  plate 
containing  an  ideal  (i.e.,  straight,  infinitely  thin,  and  perfectly  insulated)  crack.  The 
model  will  be  used  to  examine  the  effect  of  system  parameter  variability  on  flaw 
detectability.  The  results  of  this  sensitivity  analysis  will  then  be  combined  with  empirical 
noise  distributions  to  determine  probability  of  detection.  A  goal  in  this  research  is  to 
develop  a  method  where  simulation  represents  the  experimental  approach,  is  cheaper  and 
faster,  and  identifies  sources  of  unreliability  in  SQUID  NDE. 
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CHAPTER  II 


SQUID  NDE  MEASUREMENT  MODEL 

Modeling  the  Measurement  Process 

Due  to  the  large  number  of  parameters  defining  a  method,  experimental  results 
alone  are  usually  inadequate  for  examining  the  effect  of  parameter  variation  on  POD. 
Therefore,  predictive  measurement  models  are  needed  in  reliability  analysis  to  extend  the 
database  required  for  a  POD  analysis.  Using  mathematical  modeling,  the  response  of  the 
measurement  system  to  the  anomalies  of  interest  (e.g.,  cracks,  corrosion,  and  voids)  can 
be  simulated  if  enough  is  known  about  the  field-flaw  interaction  that  generates  the 
response  function.  Due  to  the  complexity  of  the  mathematics  of  these  interactions, 
idealized  models  are  normally  used  but  these  still  provide  sensitivity  analysis  information 
useful  in  constructing  POD  models. 

The  goal  of  a  measurement  model  would  be  to  predict,  through  repeated 
simulations  with  random  perturbations  in  the  system  parameters  that  have  known 
distributions,  a  distribution  of  flaw  signals  about  some  mean  value.  This  distribution 
determines  the  POD  performance  for  the  particular  system/technique  being  verified. 

Flaw  Modeling 

The  objective  of  the  SQUID  measurement  model,  which  solves  the  forward 
problem  of  calculating  the  magnetic  field  from  a  known  current  distribution,  is  to 
simulate  what  the  instrument  may  see  in  a  test  environment  on  an  vinknown  sample.  A 
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SQUID  gradiometer  measures  the  vertical  component  of  magnetic  field,  which  is  the 
result  of  currents  flowing  in  the  sample.  A  flaw  will  perturb  these  currents  and,  if  the 
perturbation  is  large  enough  (e.g.,  signal/noise  >1),  show  up  as  an  anomaly  in  the 
magnetic  field  map.  Ideally,  “inverting”  this  magnetic  field  to  give  the  underlying 
current  map  which  generated  the  field  would  give  much  detail  about  the  flaw  producing 
the  anomaly  but,  unfortunately,  inverse  solutions  are  non-unique  (in  3-D)  and  even 
obtaining  approximate  reconstructions  can  be  a  very  difficult  job.  The  specific  aging 
aircraft  problem  being  addressed  is  the  modeling  of  a  fatigue  crack  emanating  radially 
from  a  fastener  hole  (schematically  shown  in  Fig.  2.1)  as  discussed  in  the  previous 
chapter  regarding  fuselage  skin  splices.  The  measurement  model  involves  only  the 
forward  problem  and  requires  that  the  flaw  geometry  be  known  beforehand. 


Fastener  hole 


Figure  2. 1  Schematic  of  a  fastener  hole  with  crack 

Superposition  Measurement  Model 

To  model  a  hole  with  a  crack,  an  earlier  measurement  model  was  constructed 
using  existing  analytical  solutions.  Several  analytical  solutions  for  very  simple 
geometries  have  been  derived  to  represent  the  magnetic  field  produced  by  a  flaw  in  a 
conductor.  In  the  case  of  direct  current  injection,  these  include  a  spherical  flaw  in  a 
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conducting  half  space,  a  flat-bottomed  cylindrical  flaw  in  a  half  space,  and  an  elliptical 
flaw  in  a  conducting  plate  (2-D).  Scans  of  standard  flaw  specimens  have  experimentally 
validated  these  models  [1 1,12]. 

The  earlier  measurement  model  constructed  used  the  solution  of  the  perturbed 
magnetic  field  due  to  an  elliptical  flaw  in  a  conducting  sheet  containing  injected  current. 
In  the  elliptical  model,  the  major/minor  axes  can  be  varied  and  also  allows  the  limiting 
case  of  when  the  minor  axis  goes  to  zero  (approximating  a  crack).  By  using 
superposition  (Bom  approximation),  one  ellipse  can  act  as  the  fastener  hole  and  another 
ellipse,  translated  to  one  side,  can  act  as  a  crack  emanating  from  that  hole  as  shown  in 
Figure  2.2. 


Fastener  hole 


Figure  2.2  Hole/crack  superposition  model  using  ellipses 

The  superposition  is  only  an  approximation  since  it  assumes  non-interaction  between  the 
hole  and  the  crack  but  provides  a  view  of  the  critical  issues  confronting  crack  detection. 

Boundary  Element  Method  Measurement  Model 

The  Boundary  Element  Method  (BEM)  provides  the  basis  for  the  updated 
measurement  model  of  a  hole  with  a  crack  and  is  more  accurate  than  the  superposition 
model  used  for  preliminary  analysis.  The  superposition  model  was  not  accurate  for  the 
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geometry  used  since  the  hole  and  the  crack  interact  when  they  are  in  close  proximity.  The 
BEM  formulation  solves  for  the  potential  problem  of  the  hole  and  the  crack  together  and 
thus,  accounts  for  their  interaction.  The  BEM  model  has  replaced  the  superposition 
measurement  model  for  use  in  generating  signals  used  in  the  POD  analysis. 

By  applying  the  law  of  Biot  and  Savart  to  the  modeling  of  dc-current  injection,  an 
integral  for  the  magnetic  field  results  in  a  form  that  can  be  solved  using  BEM  techniques. 
The  law  of  Biot-Savart  states  that  the  magnetic  field  at  a  field  point  c  can  be  determined 
from  the  volume  integral  over  the  region  91  of  the  curl  (with  respect  to  the  point  c)  of  the 
gradient  of  the  scalar  potential  V,  at  the  source  point  q,  divided  by  the  distance  r  between 
the  source  and  field  points 


An 


r{c,q) 


dr{q)  . 


(2.1) 


But 


V,x 


r{c,q)  J  r{c,q) 


=  -V,x 


(2.2) 


and,  for  the  case  of  a  homogeneous,  isotropic  medium,  V  satisfies  Laplace’s  equation 
(i.e.,  Fis  a  scalar  potential  field  and  V  x  [v(5ca/ar)]  =  0 )  and  the  second  term  on  the  right 
hand  side  of  the  above  equation  equals  zero.  Therefore, 


’(c,q) 


/ 


:-V,X 


(2.3) 


and  now  the  integral  is  a  function  of  q  only  (so  we  can  hereafter  drop  the  q  notation). 


An  [  r 


dr. 


(2.4) 
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One  of  the  vector  integral  theorems  [13]  states  that  for  a  general  vector  F, 

\VxFdT  =  \  nxFda,  (2.5) 

JiR  Jn 

where  represents  a  2-D  surface.  This  is  derived  from  the  divergence  theorem 

f  V-Adt  =  ^  A-ftdQ  (2.6) 

hi  Jn 

by  a  substitution  of  A  =  Fxc  where  c  is  some  arbitrary  constant  vector  [14]. 

By  applying  this  vector  integral  theorem  to  Eq.  (2.4)  and  expanding  out  the  cross 
product  term,  the  integrand  term  can  now  be  represented  as 


nx- 


VV  1 


3^  3^ 

- — 

^3  ^  ^ 


3^  ^V-  (3^  3^'^ 


n. 


3 


.  (2.7) 


'‘3  ‘  3  J 

Now,  assuming  that  the  integration  over  the  volume  can  be  accomplished  by  taking  thin 
slices  perpendicular  to  the  z-axis,  then  the  resulting  surface  integral  of  this  slice  is  over  a 
top  and  bottom  surface  and  a  bounding  side  surface  (see  Fig.  2.3). 


Figure  2.3  Surface  integration  of  slice 

Looking  at  each  one  of  these  surfaces  separately: 

1)  Top  surface:  normal  is  in  positive  z-direction  only. 


«x- 


vr 


1 


-n. 


'  ^lop  '/ 

where  dV ldy,dV l3  are  evaluated  at  the  top  surface 


1  {  3^h 
r,\  3  j 


(2.8) 
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2)  Bottom  surface:  normal  is  in  negative  z-direction  only, 


«x- 


VV 


y  bottom 


n 


/+-  -n 


£^V- 

<3c  y ' 


(2.9) 


where  are  evaluated  at  the  bottom  surface.  Note  here  that  r^  and 

differ  by  the  slice  thickness  h  in  their  z-component,  i.e.  if  r,  is  represented  by 


r,  =^|{x,-  xf  +  (y,  -  yf  +  -  zf  ,  (2.10) 

where  the  subscript  c  denotes  the  coordinates  of  the  pickup  coil,  then  r^,  is 
represented  by 


(2.11) 


3)  Bounding  side  surface  -  normal  is  in  x-y  plane  only, 
VF^l  \(  If  if 


3^  ^  Ir 

(2.12) 


where  all  derivatives  are  evaluated  at  the  bounding  side  surface. 


The  total  surface  integral  is  then  the  sum  of  these  three  cases.  If  it  is  further 
assumed  that  only  2-D  current  flow  exists  within  this  slice,  then  there  will  be  no  variation 
of  the  potential  in  the  z-direction  and  the  resulting  derivatives  of  the  potential  at  the  top 
and  bottom  surfaces  will  be  equal 


3c 

3F 


=  0 


3c 

cV 


wottom 


3c 

3V 


wottom 


(2.13) 
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The  cross  product  term  is  now 


'i  'iA  ^)  I'i  nh  <^J  r^K  ^  ^) 


This  relationship  shows  that  the  magnetic  field,  for  the  assumption  of  2-D  current  flow, 
has  components  in  all  three  directions.  For  the  conventional  SQUID  magnetometer 
system,  the  pickup  coil  is  oriented  in  the  x-y  plane  and  thus  measures  only  the  z- 
component  of  magnetic  field.  Therefore,  the  only  component  of  interest  here  is 


f.  VF^ 

1 

YIX - 

=  — 

2D 

r 

\  / 

2-direction 

If  dV 

-  n, - —  \k 


(2.15) 


This  result  can  be  rewritten  as  VF/r  dotted  into  the  tangent  vector: 


[ds’ds 


,  (dy  dx]  -  . 
n=  -j--—  U  ■n  =  0 
\  ds  ds  ) 


(2.16) 


t=[-ny,n,) 

.  vv  if  3^  3^)  (vF-r) 

nx - =  -  w - - - 

r  r\.  ^  ^  dx )  r 


So,  with  this  substitution  and  letting  q  -^Q  on  the  boundary,  the  final  result  for  the 
magnetic  field  (in  the  z-direction)  at  the  field  point  c  as  a  function  of  the  tangential 
derivatives  of  the  scalar  potential  is 

Att  r{c,Q) 

For  the  assumption  of  2-D  current  flow  within  the  thin  slice  (thickness  h),  the  integrand 
of  Eq.  (2.17)  is  approximately  constant  in  the  z-direction  £ind  the  surface  integral  can  be 


represented  as  a  line  integration  multiplied  by  the  thickness  h.  In  this  way,  a  pseudo 
surface  integration  can  be  done  over  the  surface  Q  where  dQ  =  hds  (see  Figure  2.4).  The 
final  form  of  the  boundary  integral  to  be  solved  can  now  be  represented  as 


Att 


i 


r{c,Q} 


Mq)- 


(2.18) 


Figure  2.4  Relationship  between  surface  and  line  integration 


Some  care  must  be  exercised  here  to  keep  track  of  the  sign  as  you  move  around  the 
boundary  during  the  integration  of  Eq.  (2.18).  Convention  is  to  define  the  positive 
direction  along  5  as  to  keep  the  material  on  your  left  as  you  move  along  the  boundary. 
Integration  direction  around  the  boundary  of  the  fastener  hole  with  a  crack  is  shown  in 
Fig.  2.5.  For  the  upper  surface  of  the  crack,  the  integration  direction  is  in  the  positive  x 
direction  (ds  =  dx)  but,  for  the  lower  surface  of  the  crack,  the  integration  direction  in  the 
negative  x  direction  {ds  =  -dx). 
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Figure  2.5  Boundary  integration  direction 

Appendix  A  details  the  mathematical  formulation  [15]  that  will  be  applied  to 
solving  these  boundary  integral  equations.  In  this  formulation,  the  crack  is  modeled 
through  the  Green’s  function  and  not  through  the  boundary  of  the  crack  eliminating  the 
difficulties  and  inaccuracies  associated  with  mesh  construction  around  the  crack  tip 
region.  The  boundaries  of  interest  here  have  been  separated  into  those  associated  with  the 
crack  (s  =  f)  and  non-crack  (s  =S)  boundaries  since  different  techniques  will  be  required  to 
carry  out  the  integration  for  each. 

Non-crack  boundaries  (s  =  S) 

The  BEM  program  calculates  dVjdt  as  a  piecewise  linear  result  since  V(Q)  is 
piecewise  quadratic.  Letting  F,  =  VF  -f  be  the  tangential  component  of  the  gradient  of 
the  scalar  potential,  we  have  for  the  line  segment  along  s 

VXQ)  =  b,+b,s,  (2.19) 

where  b^  and  are  different  for  each  line  segment.  This  is  done  for  all  boundary  surfaces 
except  for  the  crack  surface  (T)  {e.g.,  fastener  hole,  plate  edges,  etc.).  Substituting  this  V^ 
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into  the  integral  for  B/c),  Eq.  (2.18),  gives  the  contribution  to  the  magnetic  field  at  point 
c  due  to  all  non-crack  boundaries 


ds{Q). 


(2.20) 


To  carry  out  the  integration,  the  geometrical  relations  for  r  and  s  must  be  determined. 


Figure  2.6  Geometrical  relationships  for  integration 


From  the  Fig.  2.6  we  can  see  that 

D 

r  = - 

cos^ 

,  rd9  Dd9 

dS  = - = - r - 

COS  9  cos  9 

s  =  Dtan^l^^  =>  Z)(tan0-tan^,). 
Therefore,  for  non-crack  boundaries,  the  integral  is  now 


(2.21) 
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(2.22) 


where  5,  =  6,  -  b^DimO^  (constant  term).  Since  the  surface  S  is  divided  into  n  segments, 
the  summation  of  this  equation  evaluated  over  each  segment  for  which  F,  is  defined  by  a 
single  linear  equation  (each  with  a  different  6,  and  will  determine  at  the  pickup 
coil; 


B.(c)  =  S^ 


b^D 
cos  6 


(2.23) 


where  and  6*2  correspond  to  the  angles  between  the  perpendicular  D  and  each  endpoint 
of  the  segment  as  the  summation  moves  in  a  positive  sense  around  the  boundary  (material 
on  left). 


Crack  boundaries  (s  =  1) 

For  the  crack  (length  a),  the  singular  behavior  of  V^  at  the  crack  tip  requires 
numerical  integration  to  evaluate  the  magnetic  field  contribution  of  the  crack.  Values  of 
V^  will  be  evaluated  near  the  upper  and  lower  crack  surfaces  using  the  BEM  program  and 
then  the  magnetic  field  will  be  calculated  through  a  numerical  integration  of  the 
previously  stated  Biot-Savart  relation  (Eq.  2.1).  The  series  expansion  of  F,  is  given  by 
[16] 
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Ix^-a^ 


+  /2W, 


(2.24) 


where  yj(x)and  /2(x)  are  analytic  functions.  The  first  term  contains  the  discontinuity  at 
the  crack  tip  while  the  second  term  is  needed  to  match  far  field  boundary  conditions.  On 
the  crack  surface  (/)  this  can  be  rewritten  as 

^|a- x-Ja  +  X 


_  1  ifXx) 


+  /2W 


(2.25) 


-  r—  iF\  (^)  +  yi  (^)> 

ip 

where  p  is  the  distance  measured  from  the  crack  tip  (see  Fig.  2.7).  Equation  (2.25)  is 
well  behaved  except  at  the  crack  tip.  On  the  crack  boundary  near  the  crack  tip,  the 
behavior  of  V^  can  be  represented  in  terms  of  the  potential  intensity  factor  (PIF), 

p-*0 


v\  .  =±-^+... 

^  Icrack  lip  region  //%_  _  ‘ ’ 

y2;ryO 


(2.26) 


where  the  plus/minus  sign  is  associated  with  the  upper  and  lower  crack  surfaces 
respectively.  Note  that  at  p  =  0  (crack  tip),  F,  =  00.  The  IC  factor  is  proportional  to 
crack  length  (square  root)  and  also  to  current  and  can  be  thought  of  as  representing  how 
much  current  is  being  diverted  around  the  tip  of  the  crack;  the  longer  the  crack  length,  the 
more  current  that  is  diverted.  By  combining  Eqs.  (2.25)  and  (2.26),  we  can  express  F,  at 
the  crack  tip  in  terms  of  the  PIF 


(2.27) 


21 


Figure  2.7  Distance  p  measured  from  crack  tip 

It  will  be  assumed  that  /2(a),  the  far  field  boundary  term,  is  small  relative  to  the  crack  tip 
term  and  can  be  neglected.  With  the  function  value  at  the  crack  tip  now  defined,  we  can 
continue  with  the  numerical  integration  development. 

Numerical  integration 

Numerical  integration  over  the  crack  is  accomplished  through  discretization  of  the  crack 
boundary  into  elements.  By  utilizing  a  coordinate  transformation,  F,  can  be  expressed  in 
terms  of  nodal  values  and  interpolation  functions  (shape  functions)  of  an  intrinsic 
coordinate  Once  mapped  into  ^-space,  Gaussian  quadrature  numerical  integration  can 
be  used  to  evaluate  the  integral.  This  formulation  uses  quadratic  interpolation  functions 
Nm(^)  which  requires  three  nodes  per  boundary  element.  The  crack  is  modeled  with  one 
crack-tip  element  containing  the  crack-tip  node  and  one  regular  element  (Fig.  2.8). 

The  regular  element  is  a  straightforward  application  of  the  interpolation  functions 
but  the  crack-tip  element  will  use  a  quarter-point  formulation  [17]  to  map  the  singular 
behavior  at  the  crack  tip.  The  general  forms  of  the  quadratic  interpolation  functions  used 
to  map  the  function  to  be  integrated  into  the  ^-space  are 
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crack-tip  element 


S5  s. 


S, 


regular  element 


s,  s, 
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Figure  2.8  Discretization  of  crack  into  elements 


NSt)=\4(4-\) 

(2.28) 

A',fe)=iffe+1), 

which  are  shown  in  Fig.  2.9.  The  coordinate  transformation  into  ^-space  of  the  integral 
results  in  the  following  sum 


fF(j)&  ,  (2.29) 

r  *=i  -I 

where  is  the  value  of  the  function  at  the  node  and  J  is  the  Jacobian  of  the  mapping 


that  accounts  for  the  spatial  scaling  associated  with  the  mapping 


ds  dN^ 

d^  h  d^  *■ 


(2.30) 


Figure  2.9  Quadratic  interpolation  functions  (a)  regular  element  (b)  crack-tip  element 
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For  those  boundary  elements  (length  L)  that  do  not  contain  the  crack  tip,  the  resulting 
integral  for  becomes 
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(2.31) 


which  can  now  be  numerically  integrated  using  Gaussian  quadrature  techniques. 

For  the  crack-tip  element  on  the  upper  crack  surface,  the  standard  mapping  to 

space  needs  to  be  modified  to  accommodate  the  l/ ^fs  behavior  of  V,(s)  at  the  crack  tip 
(see  Fig.  2.9  (b)).  The  general  approach  is  to  split  the  1^(5')  term  into  a  singular  term 
multiplied  by  a  non-singular  coefficient.  The  form  of  V,  (s)  follows  from  Eq.  (2.25) 

if{s)  1 


(2.32) 


To  map  V,{s)  into  ^-space,  the  singular  term  uses  an  inverse  mapping  relationship  and 
the  non-singular  term  uses  the  quadratic  shape  functions.  For  the  inverse  mapping, 

the  l/ yfs  behavior  is  represented  by  placing  the  midnode  of  the  quadratic  element  at  the 
quarter  point  location 


#  =  -1  +  2 


fi 


^fs  (# + 1) 


(2.33) 


The  Jacobian  is  given  by 
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(2.34) 


Note  that  at  the  point  where  =  0  (^  =  -1),  the  Jacobian  also  equals  zero.  This 


characteristic  is  important  when  evaluating  the  integral  of  the  mapped  function. 
Therefore,  the  5^  contribution  due  to  the  crack-tip  element  on  the  upper  crack  surface  is 


(2.35) 

f  (#  +  !)<(# 


Also  note  that  as  the  singularity  term  (4+1  in  denominator)  goes  to  infinity  but  is 

canceled  by  the  Jacobian  term  (4+1  in  numerator)  going  to  zero,  thereby  making  the 
overall  function  finite  in  the  mapped  space.  All  values  of  F,  r,  and  5  are  known, 

including  F(s^)  =  K* I -Jin  (from  Eq.  2.27)  at  the  crack  tip  node  that  has  the  singularity. ' 
The  crack-tip  numerical  integration  using  the  quadratic  interpolation  functions  and  the 
quarter-point  formulation  used  in  the  model  was  verified  against  an  integrable  test  case 
having  the  same  singular  behavior  as  Eq.  (2.24) 
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Using  the  Measurement  Model 


A  computer  program  (Appendix  B)  has  incorporated  this  BEM  development  into 
a  simulated  scan  of  a  SQUID  over  a  dc-current  injected  plate  containing  a  crack.  The 
program  can  accommodate  various  sample  geometries  and  can  include  an  interior 
boundary  (e.g.,  a  fastener  hole)  if  desired.  Figure  2.10  shows  the  basic  input  variables. 

Calculation  of  the  magnetic  field  is  a  two  step  process.  First  the  BEM 
subprogram  must  be  run  to  solve  for  the  potentials  needed  for  the  integrals  involved  in 
the  magnetic  field  calculation.  This  step  requires  that  a  boundary  mesh  be  generated 
reflecting  the  geometry  of  the  sample.  The  crack  is  handled  through  the  Green’s  function 
formulation  and  does  not  need  to  be  modeled  in  this  mesh. 


User 
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gradiometer  coH  radii 


gradiometer  Hftoff 
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Mode! 


Figure  2.10  Input  variables  for  BEM  measurement  model 
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Once  the  BEM  subprogram  has  calculated  the  required  nodal  tangential  derivatives,  the 
measurement  model  then  steps  through  each  x-y  scan  coordinate  and  calculates  the 
magnetic  field  contributions  from  the  tangential  derivatives  at  all  the  boundary  nodes 
(crack  and  non-crack).  The  magnetic  field  is  integrated  over  the  areas  of  the  gradiometer 
coils  to  determine  net  flux.  This  output  is  then  plotted  in  the  form  of  a  2-D  magnetic  field 
map. 

If  the  measurement  model  is  to  be  used  to  generate  POD  analysis  data,  the 
program  is  set  up  to  loop,  using  the  same  BEM  input  file,  while  sampling  from  a  noise 
(or  system  parameter)  distribution  in  Monte  Carlo  fashion.  The  peak-to-peak  value  of  the 
magnetic  field  is  extracted  from  each  2-D  scan  (one  Monte  Carlo  loop)  and  represents  the 
signal  used  in  the  determination  of  POD.  In  this  way,  a  distribution  of  signals  can  be 
generated  for  a  single  crack  length.  For  each  crack  length,  the  BEM  subprogram  is  run 
only  once,  then  the  signal  distribution  for  that  crack  generated  from  the  measurement 
model  Monte  Carlo  loops.  Calculation  of  the  POD  from  this  data  is  discussed  in  Chapter 
IV. 
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CHAPTER  III 


COMPARISON  OF  BEM  MEASUREMENT  MODEL  TO  EXPERIMENT 

Comparison  Issues 

Several  experimental  measurements  were  compared  to  the  BEM  measurement 
model.  There  are  some  important  issues  affecting  how  closely  the  measurement  model 
can  represent  the  experimental  measurement.  Fabrication  of  test  samples  to  be  used 
specifically  for  comparison  to  models  is  usually  a  compromise  between  representing  the 
ideal  geometries  used  in  the  calculations  and  what  is  found  in  an  actual  aircraft  structure. 
In  addition,  the  uncertainties  associated  with  taking  experimental  measurements  make 
comparison  as  much  a  qualitative  analysis  as  it  is  a  quantitative  one.  These  issues,  and 
others,  will  be  discussed  in  more  detail  in  Chapter  V. 

Also,  direct  amplitude  comparison  requires  accurate  SQUID  calibration  factors. 
A  description  of  a  calibration  procedure  used  in  this  laboratory  is  given  in  Appendix  C. 
Only  one  of  the  following  experimental  comparisons  utilized  this  factor. 

Experimental  Measurements 

Experimental  measurements  were  made  using  three  different  SQUID  systems  for 
comparison  to  the  BEM  measurement  model.  These  experimental  measurements  do  not 
necessarily  validate  the  BEM  measurement  model  but  show  that  the  model  does  do  well 
in  determining  the  magnetic  fields  associated  with  common  test  conditions  found  in  the 
laboratory. 
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Calibrated  System:  SBIR 

A  new  optimization  procedure  was  conducted  for  one  channel  of  a  5-channel  low- 
Tc  SQUID  gradiometer  system  (termed  SBIR  SQUID)  to  determine  an  accurate  value  of 
the  voltage-to-magnetic-field  factor  required  to  convert  the  output  voltage  of  the  SQUID 
to  units  of  tesla,  which  is  consistent  with  the  measurement  model.  This  was  done  so  that 
at  least  one  comparison  based  on  accurate  experimental  signal  amplitudes  could  be  done. 
The  other  SQUIDs  used  in  these  measurements  do  not  have  optimized  calibration  factors 
at  this  time. 

SBIR  SQUID  is  documented  as  having  an  axial  gradiometer  with  5-mm  diameter 
pickup  and  balance  coils  with  a  baseline  of  20  mm.  The  calibration  factor  that  was 
determined  for  channel  1  of  SBIR  SQUID  is  assumed  to  be  independent  of  magnetic  field 
source  making  it  applicable  to  scanning  samples  with  different  associated  magnetic 
fields.  A  measurement  was  made  using  a  100  mA  dc  current  injected  into  a  100  mm  x 
150  mm  X  0.03  mm  copper  clad  circuit  board  containing  a  9  mm  diameter  hole  with  a  9 
mm  X  0.03  mm  slot  on  one  side  of  the  hole  (see  Figure  3.1).  Although  this  setup  does  not 
provide  completely  uniform  current  injection  across  the  sample  (transverse  to  the  slot) 
due  to  the  point  source  electrodes,  the  region  around  the  central  region  should  be 
relatively  uniform. 
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100X150X0.03  mm  Cu  plate 


Figure  3.2  shows  the  contour  map  resulting  from  a  2-D  scan  (unshielded 
environment)  over  the  sample.  The  optimized  parameters  (pick-up  coil  radius  =  2.8  mm, 
baseline  =  2.06  cm,  and  liftoff  =  14  mm)  were  used  in  the  BEM  measurement  model  and 
both  magnetic  maps  compared.  Figure  3.2  also  shows  the  profile  comparison,  after 
scaling  the  experimental  result  with  the  optimized  calibration  factor  (6.16  x  10‘®  TfV), 
corresponding  to  the  line  AA'  on  the  contour  map.  At  this  large  of  a  liftoff,  the  dipole 
signal  does  not  reveal  much  about  the  hole/slot  geometry  but  the  interest  here  was  to 
ensure  that  the  measurement  model  was  calculating  field  amplitudes  correctly.  Using  the 
optimized  parameters,  the  measurement  model  is  within  2%  of  experiment  (measuring 
Bpp  of  the  dipole  associated  with  the  hole  and  slot).  Some  difference  can  be  attributed  to 
the  assumption  of  uniform  current  density  across  the  plate,  which  can  have  a  large  effect, 
and  is  most  likely  causing  the  amplitude  error  seen  at  both  ends  of  the  scan  profile.  The 
current  density  non-uniformity  can  also  be  seen  on  the  experimental  contour  map  as 
asymmetries  in  the  contour  lines.  Errors  in  the  optimized  parameters  or  calibration 
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number  would  show  up  as  a  mismatch  of  the  profiles  in  shape  and  amplitude  beyond  that 
attributable  to  experimental  error. 


Figure  3.2  SBIR  magnetic  field  contour  map  and  comparison  profile 


Measurements  Using  Other  Systems 

Other  comparisons  were  made  using  SQUID  systems  that  did  not  have  optimized 
calibration  factors.  By  scaling  the  amplitudes  of  the  measurement  model  map  to  the 
experimental  map,  the  use  of  the  calibration  factor  was  bypassed  and  then  only  shape 
comparisons  were  made. 
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MicroSQUID 

Using  MicroSQUID,  a  comparison  experiment  used  a  5  mA  dc  current  injected 
into  a  75  mm  x  150  mm  x  0.03  mm  copper  clad  circuit  board  containing  a  15  mm  x  0.03 
mm  slot  cut  with  a  scalpel  (see  Figure  3.3).  Figure  3.4  shows  the  contour  map  resulting 


75X150X0.03  mm  Cu  plate 


from  a  2-D  scan  (shielded  environment)  over  the  sample  with  the  lines  AA',  BB',  and 
CC'.  Profiles  along  these  lines  on  this  map  are  compared  to  those  calculated  by  the 
measurement  model.  Magnetic  field  shape  characteristics  (matching  of  peaks  andvalleys) 
reflect  the  accuracy  of  the  measurement  model  since  a  difference  between  the  model  and 
experiment  would  show  up  as  a  mismatch  of  the  profiles  at  either  the  edge  or  the  crack 
locations.  The  measurement  model  is  in  very  good  agreement  (only  shape  characteristics 
compared)  with  experiment  for  all  profiles.  A  slight  mismatch  at  the  right  side  of  AA' 
can  be  seen  and  is  most  likely  due  to  a  small  variation  in  liftoff  during  a  scan  (i.e.  sample 
not  level). 
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Figure  3.4  MicroSQUID  magnetic  field  contour  map  and  comparison  profiles 


Mobile  HTS  SQUID 

Another  measurement  used  a  50  mA  dc-current  injected  into  a  100  mm  x  150  mm 
X  0.03  mm  copper  clad  circuit  board  containing  a  9  mm  diameter  hole  with  a  9  mm  x 
0.03  mm  slot  on  one  side  of  the  hole  (see  Figure  3.5).  However,  the  SQUID  used  for  this 
measurement  was  a  high-Te  system  [24]  which  used  a  planar  gradiometer  (the  pickup  and 
balance  coils  are  in  a  plane  parallel  to  the  sample  surface)  instead  of  ein  axial  one.  Figure 
3.6  shows  the  contour  map  resulting  from  a  2-D  scan  (imshielded  environment)  over  the 
sample.  The  quadrapolar  shape  results  from  the  planar  gradiometer  acting  as  a  spatial 
differentiater  and,  by  taking  the  derivative  of  the  dipole-shaped  magnetic  signal,  results  in 
a  quadrapole  shape.  Since  the  planar  gradiometer  was  oriented  parallel  to  the  plate 
edges,  the  edge  signal  is  approximately  zero  and  therefore,  only  that  part  of  the  scan 
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1 00  X 1 50  X  0.03  mm  Cu  plate 


above  the  hole/slotwas  used  in  these  comparisons.  Figure  3.6  also  shows  the  profile 
comparisons,  after  scaling,  corresponding  to  the  lines  AA',  BB',  and  CC'  on  the  contour 
map.  The  measurement  model  is  in  very  good  agreement  (only  shape  characteristics 
compared)  with  experiment  for  all  profiles.  Again,  a  difference  between  the  model  and 
experiment  would  show  up  as  a  mismatch  of  the  profiles  at  the  location  of  the  hole/crack. 


A 

— ^ — . — I — I — 1 — 1 — 1 — I — 1 — 1 — ■ — I— 

-  A-A' 

— _ _ I  . . . . . L- _ _ _ _J 

Figure  3.6  Mobile  HTS  magnetic  field  contour  map  and  comparison  profiles 


34 


CHAPTER  IV 


POD  CURVES  FROM  SIMULATED  EXPERIMENTAL  DATA 


Probability  of  Detection 


Introduction 

In  general,  inspection  systems  must  use  some  criterion  to  “accept  or  reject”  the 
part  being  tested.  A  threshold  (associated  with  a  minimum  crack  size)  can  be  set  so  that 
the  signals  above  the  threshold  are  rejected  and  those  below  accepted.  Figure  4.1  shows 
the  cumulative  distribution  of  rejected  signals  for  a  range  of  signals  around  the  threshold. 
For  an  ideal  inspection  system,  no  signals  below  the  threshold  would  be  rejected  and 
100%  of  those  above  the  threshold  would  be  rejected.  But,  for  real  systems,  uncertainties 
in  the  measurement  process  cause  the  curve  to  be  rounded,  creating  two  regions  of  error. 
False  accepts  (type  I  error)  are  signals  above  the  threshold  that  are  missed  leading  to 
issues  of  safety.  False  rejects  (type  II  error)  are  signals  below  the  threshold  that  are  being 
unnecessarily  rejected  leading  to  issues  of  cost  of  repair  or  replacement. 


Figure  4.1  Cumulative  distribution  of  rejected  signals:  real  versus  ideal 
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Determination  from  Experimental  Data 

Determination  of  POD  requires  that  the  underlying  signal  distribution  of  the  data  be 
known  (or  assumed).  The  distribution  of  signals  results  from  such  things  as  sample  related 
variability  (e.g. ,  material  properties,  crack  geometry  and  orientation)  and  measurement 
system  variability  (e.g.,  sensor  characteristics,  liftoff,  scanning,  calibration,  and  different 
users).  The  POD(a)  is  determined  by  integrating  this  signal  distribution,  fa{a),  above 
some  specified  signal  threshold  a,^ , 

00 

POD(a)=  |f,(a)dfl  (4.1) 

au. 

as  represented  by  the  shaded  areas  in  the  upper  part  of  Fig.  4.2.  This  integration  is  done  at 
each  crack  length  resulting  in  the  POD  curve  shown  in  the  lower  part  of  Fig.  4.2. 


Figure  4.2  Integration  of  signal  to  determine  POD 
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If  the  analysis  is  vising  simulated  experimental  data  through  use  of  the  BEM 
measurement  model  and  Monte  Carlo  methods,  then  the  POD  at  a  given  erack  length  is 
determined  through  direet  numerical  summation  of  those  signals  above  the  threshold 

POD{a)  oc  ^  (response  signals  >a,^).  (4 .2) 

N 

This  approach  does  not  require  that  the  distribution  of  the  signals  be  known. 

For  real  experimental  data,  it  is  commonly  assumed  that  the  distribution  of  signals 
at  a  fixed  crack  length,  /„(«),  is  normally  distributed.  Then  a  functional  form  is  used  to 
determine  POD.  Previous  studies  [25]  have  shown  that  NDI  experimental  data  commonly 
have  a  POD  functional  form  that  follows  the  log-odds  or  log  logistic  model 


PODia)  = 


1  +  exp 


7rf  Ina-//' 


(4.3) 


where  //represents  the  central  location  of  the  POD  curve  and  cris  a  measure  of  the  flamess 
of  the  POD  curve  (larger  a,  flatter  POD  curve).  Another  POD  function  that  has  been  used 
[26]  is  the  Weibull  function 

POD(a)  =0  a<a,^ 


=  1  -  exp 


a 

a 

[  V  J 

a>a 


th  9 


(4.4) 


where  a  is  the  bandwidth,  and  rj  is  the  central  tendency  of  the  POD  function. 

NDE  reliability  data  normally  falls  into  two  types,  either  data  in  which  the 
inspection  result  was  recorded  as  pass  or  fail,  meaning  that  a  flaw  was  either  found  or  not 
found,  or  data  in  which  the  measuring  instnunent  signal  d,  in  response  to  an  actual  crack  of 
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size  a,  was  recorded.  For  pass/fail  data,  the  distributional  parameters  {e.g.,  //  and  cr)  for  the 
POD  function  can  be  determined  using  maximum  likelihood  methods.  Details  describing 
these  methods  can  be  found  in  [25, 27].  For  signal  data,  regression  analysis  can  be  used  to 
estimate  the  distributional  parameters.  In  this  case,  the  log-odds  model  is  reparameterized 
to 


A 

which  shows  that  the  data  plots  linearly  in  log-space  and  a  and  /?  can  be  determined  from 
linear  regression. 

SQUID  systems  measure  signal  type  data  as  an  output  voltage  proportional  to  the 
magnetic  field  being  measured.  Signal  data  has  the  advantage  over  pass/fail  data  since  it 
contains  more  information  and,  as  a  result,  requires  a  smaller  minimum  number  of  samples 
to  conduct  an  analysis.  For  example,  it  is  recommended  [25]  that  at  least  60  flaw  samples 
be  used  for  a  pass/fail  analysis  versus  30  for  an  d  analysis.  Other  advantages  of  signal  data 
include: 

-  The  underlying  assumptions  of  the  POD  model  are  easily  tested  with  statistical  tools. 
Since  the  POD  model  is  derived  from  the  correlation  of  the  d  versus  a  data,  the 
assumptions  concerning  the  POD(a)  model  can  be  tested  using  the  signal  data.  Also, 
the  pattern  of  d  shows  the  aeceptable  range  of  extrapolation. 
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-  The  decision  threshold  is  set  after  the  signals  are  recorded  allowing  examination  of  the 
effect  of  different  settings  on  POD  without  rerunning  the  experiment, 

-  The  POD  model  is  valid  over  large  flaw  size  range  and  it  is  not  critical  to  obtain  a  larger 
sampling  of  flaws  around  the  area  of  steep  slope  on  the  POD  curve. 


Confidence  Bounds 

Aircraft  requirements  specify  that  a  95%  confidence  bound  region  be  determined 
for  these  curves.  For  experimental  data  using  the  log-odds  model,  a  method  from  [29]  can 
be  adopted  to  calculate  the  lower  confidence  bound.  The  y-percent  lower  confidence  bound 
is  given  by 

P6>D(a)  =  (D(zJ,  (4.7) 

where 


and 


Zi  =  Z-, 


N 


ssx 


X  In  a 

—  1  ^ 

Z  =  — Tina, 

N  =  sample  size 

/I  =  y percentile  of  a  -  distribution  (2  dof) 
SSX 

/=]  JV  J 
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For  simulated  data  using  the  BEM  measurement  model,  Monte  Carlo  simulations 
will  be  used  to  determine  the  95%  lower  confidence  bound  and  only  require  that  the 
signal  and  noise  distributions  be  known.  This  will  be  discussed  further  when  these  values 
are  calculated  in  a  later  section. 

Importance  of  Setting  Threshold  Value 

It  is  important  to  note  the  dependence  of  POD  on  threshold  value.  The  threshold 
value  is  primarily  a  function  of  background  noise  and  the  choice  of  the  minimum  signal- 
to-noise  ratio  for  which  a  crack  can  still  be  seen.  This  is  very  dependent  on  the 
instrument,  technique,  and  user  of  an  NDE  inspection  system.  The  tradeoff  with  setting 
the  threshold  lower  to  increase  POD  for  a  particular  crack  length  is  the  simultaneous 
increase  of  the  probability  of  false  alarms  (POFA)  arising  from  the  background  signal 
being  mistaken  for  a  crack  signal.  The  POFA  is  calculated  in  the  same  manner  as  POD 

with  the  noise  distribution  being  integrated  above  the  threshold  value. 

00 

POFA  =  jNoise(a)  da  (4.8) 

aih 

Figure  4.3  shows  schematically  the  normalized  distributions  of  the  background 
noise  and  the  crack  signal  pointing  out  the  regions  of  concern  and  shows  the  tradeoff 
between  POD  and  POFA  when  decreasing  the  threshold. 
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Figure  4.3  Effect  of  setting  threshold  value 


Progression  of  POD  Approaches 

The  following  discussion  of  some  of  the  earlier  POD  approaches  leads  up  to  the 
present  approach  and  points  out  some  of  the  critical  items  that  were  important  in  deciding 
how  best  to  determine  POD  for  SQUID  systems. 

Failure  of  Superposition 

Initial  POD  work  utilized  a  superposition  model  consisting  of  an  elliptical  crack 
superimposed  on  a  hole  in  an  infinite  plate  (refer  back  to  Fig.  2.2).  Superposition  implies 
that  the  hole  and  the  crack  are  behaving  independently  and  do  not  influence  each  other 
when  in  close  proximity.  The  errors  associated  with  superpostion  approximations  can  be 
illustrated  by  direct  comparison  to  the  BEM  model.  Figure  4.4  displays  centerline 
profiles  for  a  5  mm  hole  with  a  5  mm  crack  at  liftoffs  of  3  mm  and  1  mm  revealing  how 
the  superposition  approximation  breaks  down  as  liftoff  distance  is  reduced. 
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Figure  4.4  Comparison  of  BEM  and  superposition  models 

It  is  apparent  that  in  these  cases  the  hole-crack  interaction  cannot  be  neglected.  However, 
for  greater  liftoffs  and/or  smaller  crack  lengths,  the  superposition  model  error  is  not  as 
significant  and  may  be  used  in  calculating  approximate  fields.  Figure  4.5  shows  that  for  a 
smaller  crack  size  of  3  mm,  the  superposition  profile  is  not  as  distorted  as  compared  to 
the  5-mm  crack  length  case  in  the  upper  part  of  Fig.  4.4. 
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Figure  4.5  Comparison  for  a  5-mm  hole  with  a  3-mm  crack 
Standard  Hole  Subtraction  Technique 

One  analysis  technique  that  was  examined  was  that  of  subtracting  a  standardized 
hole  from  the  hole-with-crack  signal  to  extract  that  part  of  the  signal  due  to  the  crack  (see 
Figure  4.6).  This  was  done  by  subtracting  the  2-D  magnetic  map  associated  with  a  scan 
over  a  standard  fastener  hole  from  the  2-D  map  of  the  hole  with  crack.  The  resulting 
“residual”  2-D  map  could  then  be  thought  of  as  the  magnetic  field  due  to  the  crack  alone. 
The  peak-to-peak  value  of  this  residual  map  was  defined  as  the  crack  signal. 


Hole  +  Crack  Hole  Crack 

Figure  4.6  Subtraction  of  standard  hole  to  extract  signal  due  to  crack 
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Unfortunately,  in  real  experimental  data,  it  is  difficult  to  determine  exactly  where 
in  the  resulting  magnetic  map  the  hole  is  located  so  that  there  will  be  an  uncertainty 
associated  with  this  subtraction  process.  Figure  4.7  shows  the  sensitivity  of  the  crack 
signal  to  the  variation  of  the  position  of  the  subtraction  hole.  The  left  graph  represents  the 
crack  signal  (the  “residual”  signal)  as  the  subtracted  hole  position  was  varied  in  a 
direction  transverse  to  the  crack.  The  right  graph  is  determined  the  same  way  with  the 
subtraction  hole  position  varying  this  time  parallel  to  the  crack.  The  relative  magnitudes 
of  the  crack  signals  for  the  two  graphs  show  that  the  response  surface  was  more  sensitive 
to  variation  in  the  subtraction  hole  position  in  the  direction  parallel  to  the  crack. 


Figure  4.7  Crack  signal  dependence  on  subtraction  hole  position  (Note  the  differing 
vertical  scales) 


Since  the  residual  crack  signal  was  dependent  on  the  position  of  the  subtraction  hole,  the 


POD  was  dependent  upon  the  uncertainty  associated  with  this  subtraction  process.  For 


illustrative  purposes,  a  distribution  was  placed  on  the  subtraction  hole  position  to 


simulate  the  experimental  situation.  An  independent  bivariate  normal  distribution. 
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A^(0, 1  mm),  was  assumed  for  the  subtraction  hole  position,  then  a  sampling  of  1000 
random  positions  results  in  a  distribution  of  residual  crack  signals  shown  histogram 
format  in  Figure  4.8. 


Figure  4.8  Histogram  showing  distribution  of  residual  crack  signals 

The  characteristics  of  this  distribution  arose  from  the  behavior  of  the  residual 
crack  signal  (variation  shown  in  Figure  4.7)  and  the  assumed  distribution  of  the 
subtracted  hole  position.  The  sharp  truncation  on  the  left  side  of  the  distribution  was  a 
result  of  the  residual  crack  signal  having  a  local  minimum  around  the  actual  hole  center; 
any  movement  away  from  this  center  resulted  in  an  increase  in  the  residual  crack  signal. 
The  peak/tail  characteristic  can  be  attributed  to  the  assumed  normal  distribution  of  the 
subtraction  hole  position.  If  the  assumed  distribution  had  a  larger  than  1  mm  standard 
deviation,  the  residual  crack  signal  distribution  would  spread  out  more  to  the  right  but  the 
truncation  would  remain  fixed. 
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The  POD  at  each  crack  length  was  determined  by  integration  of  the  residual  crack 
distribution,  determined  for  each  crack  length,  above  a  threshold  as  previously  stated  in 
Eq.  (4.1).  The  threshold  was  set  to  satisfy  SNR>2. 

The  POD  curve,  a  result  of  direct  numerical  summation  of  the  Monte  Carlo  values, 
is  shown  in  Fig.  4.9.  It  is  common  to  assume  that  the  crack  signal  distribution  is  Gaussian 
(normal).  As  a  comparison,  the  distributions  (like  the  one  shown  in  Fig.  4.8)  were  fit  with 
Gaussian  distributions  and  a  second  POD  curve  determined,  (also  shown  in  Fig.  4.9).  The 
correct  POD  curve  (shown  darker  on  the  graph)  has  a  sharp,  asymmetrical  shape.  The  sharp 
transition  to  a  POD=l  comes  from  the  narrowness  of  the  crack  signal  distribution  and  the 
asymmetry  comes  from  the  long  narrow  tail  of  this  distribution.  As  can  be  seen,  the  normal 
distributions  commonly  assumed  for  these  calculations  might  lead  to  significant  error  when 
calculating  POD.  The  crack  length  values  are  retained  for  scale  in  this  comparison  and 
should  not  be  used  to  infer  SQUID  performance. 


Figure  4.9  Calculated  POD  curves  for  hole  subtraction  scheme 
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Present  POD  Approach 

Since  the  hole  subtraction  process  added  additional  uncertainty  that  negatively 
affected  POD,  it  was  desirable  to  be  able  to  determine  POD  without  this  subtraction.  It 
has  been  determined  that  POD  can  be  calculated  from  the  hole-with-crack  signal  directly. 
For  this  approach,  the  “signal”  is  no  longer  the  peak-to-peak  residual  signal  after 
subtraction,  but  the  overall  hole  with  crack  peak-to-peak  signal.  Figure  4.10  shows  the 
centerline  profiles  of  the  magnetic  maps  resulting  from  simulated  scans  over  a  hole  with 
crack  for  four  different  crack  lengths  while  keeping  the  hole  diameter  constant. 


Figure  4. 10  Scan  centerline  profiles  for  various  crack  lengths  (with  hole) 

As  can  be  seen,  for  increasing  crack  length,  the  signal  amplitude  increases.  Also,  the 
asymmetry  increases,  pulling  the  crack  side  of  the  dipolar  signal  to  the  right  while  the 
hole  side  remains  relatively  fixed.  The  important  characteristic  for  POD  is  the  increasing 
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peak-to-peak  value  as  a  function  of  crack  length,  which  is  the  “signal”  used  in  the  POD 
analysis.  The  effect  of  increasing  crack  length  on  signal  is  apparent  without  subtracting  a 
standard  hole.  For  comparison,  Fig.  4.1 1  shows  the  corresponding  scan  profiles  for 
cracks  alone  with  no  fastener  hole  (normalized  units  of  signal  are  on  the  same  scale). 
Here,  the  profiles  are  symmetrical  but  still  show  the  important  characteristic  of  increasing 
amplitude  for  increasing  crack  length.  Either  type  of  data  (Fig.  4.10  or  4.1 1)  can  be  used 
to  calculate  POD, 
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Figure  4.1 1  Scan  centerline  profiles  for  various  crack  lengths  (no  hole) 

When  the  peak-to-peak  signals  from  Figs.  4.10  and  4.1 1  are  plotted  as  a  function  of  crack 
length,  comparisons  can  be  made  between  the  two  (see  Fig.  4.12)  at  a  given  liftoff  (each 
liftoff  distance  would  have  a  corresponding  curve).  The  crack-only  case  can  be  shifted  by 
adding  the  signal  from  a  hole  with  no  crack  (corresponding  to  the  hole-with-crack  case  at 
crack  length  =  0;  in  this  case,  a  signal  value  w  1 .5), 
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Figure  4.12  Signals  as  a  function  of  crack  length 

Looking  at  this  shifted  curve,  it  is  important  to  notice  that  the  addition  of  a 
“standard”  hole  to  the  crack  signal  is  not  equivalent  to  the  corresponding  hole  with  crack. 
The  unusual  characteristic  of  Fig.  4.12  is  that  the  BEM  model  results  in  signals  that 
increase  more  quickly  as  a  function  of  crack  length  than  with  the  cracks  alone.  This  is 
saying  that  it  is  easier  to  see  a  3  mm  crack  on  a  fastener  hole  than  it  is  to  see  a  3  mm 
crack  by  itself.  This  also  points  out  that  the  hole  subtraction  scheme  discussed  earlier 
will  not  result  in  signals  representative  of  cracks  alone  and  would  only  work  if  the 
principle  of  superposition  applied  where  the  hole  and  the  crack  behave  independently. 
Instead,  the  hole  and  the  crack  are  not  independent,  in  fact,  the  existence  of  the  hole 
seems  to  enhance  the  detection  of  the  crack.  It  has  been  suggested  that  the  hole  may  be  a 
“crack  amplifier”  due  to  the  hole’s  influence  on  the  current  density  in  the  region  near  the 
crack  causing  the  larger  signal  as  compared  to  the  signal  from  a  crack  alone  of  the  same 
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size.  This  can  be  illustrated  by  looking  at  the  current  flow  (potential  streamlines)  using 
the  relation  (normalized  to  the  free  stream  current  density) 


f  ^ 

1-4 

V  y 


smuj . 


(4.9) 


These  were  plotted  (constant  values  of  ^y^)  for  a  region  around  the  fastener  hole  in  Fig. 
4.13  and  shows  the  current  density  concentration.  For  current  injection  transverse  to  the 
crack,  the  crack  is  located  directly  in  the  region  of  high  current  density  possibly  causing 
the  greater  signal 


Figure  4.13  Current  flow  concentration  around  fastener  hole 


contribution  observed  in  the  hole-with-crack  case.  The  variation  in  current  density  can  be 
determined  by  differentiating  Eq.  (4.9)  in  this  region  {x  =  0)  with  respect  to  the  radial 
direction  transverse  to  the  current  flow. 
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For>'  >  R ,  Fig.  4.14  shows  shows  that  the  current  density  doubles  at  the  surface  of  the 
hole  and  approaches  one  away  from  the  hole. 


Figure  4. 14  Current  density  concentration  near  fastener  hole  (R  =  2.5  mm) 


One  item  to  mention  here  is  how  this  may  affect  the  superposition  model  used 
previously.  Since  the  current  density  concentration  was  not  taken  into  account  in  that 
model,  the  contribution  of  the  crack  (thin  ellipse)  may  have  been  underestimated.  If  we 
refer  back  to  Fig.  4.5,  we  can  see  that  the  right  side  of  the  superposition  curve 
undershoots  the  BEM  result.  By  considering  the  current  density  concentration,  this  error 
may  be  less.  But,  it  is  not  simply  a  multiplicative  factor  associated  with  Fig.  4.14  (e.g.,  a 
factor  of  two  near  the  fastener  hole)  since  this  would  scale  the  entire  dipole  signal 
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associated  with  the  crack  and  the  superposition  would  then  overshoot  the  BEM  result  on 
the  left  side  of  the  profile  curve.  It  might  be  possible  to  determine  some  kind  of 
weighting  function  along  the  length  of  the  superposition  crack  (heavier  at  the  end  farthest 
from  hole  down  to  a  value  of  1  at  the  end  closest)  but,  this  is  simulating  something  more 
like  a  “crack-tip  current  concentration  factor”  and  this  is  already  handled  with  the  BEM 
formulation.  However,  approximations  could  be  made  for  a  limited  range  of  uses  if  so 
desired. 

For  illustrative  purposes,  a  representative  signal  distribution  and  noise  threshold 
were  selected  and  applied  to  both  the  hole-with-crack  and  crack-alone  cases.  Figure  4.15 
shows  Fig.  4.12  with  the  appropriate  thresholds  and  signal  distributions  sketched  in.  In 
this  way  the  effect  of  the  current  density  concentration  on  POD  can  be  examined  since  it 
is  assumed  that  both  cases  have  the  same  signal  and  noise  distributions.  It  can  be  seen 
that  the  standard  hole  signal  level  defines  the  new  baseline  for  the  hole-with-crack  case. 
Experimental  techniques  exist  [30]  which  can  measure  both  the  maximum  signal  from  the 
hole  with  crack  and  the  baseline  signal  simultaneously.  However,  in  the  case  of  a  real 
crack,  current  injected  parallel  to  the  crack  will  still  result  in  a  signal  contribution  from 
the  crack  (versus  zero  contribution  in  the  BEM  measurement  model)  so,  the  baseline 
signal  from  a  hole  alone  will  have  an  associated  uncertainty  affecting  probability  of 
detection.  The  SNR  threshold  is  the  same  for  both  of  the  cases  since  the  SNR  is  a 
function  of  the  measurement  instrument,  not  the  sample.  Also,  the  signal  distribution  at  a 
fixed  crack  length  would  be  the  same  for  both  cases  unless  the  distribution  was  dependent 
on  the  fastener  hole  {e.g.,  uncertainty  in  size  or  roundness).  Usually,  this  signal 
distribution  is  dependent  on  uncertainties  in  parameters  common  to  both  samples 
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Figure  4.15  Determining  POD  from  signal  curves 

(e.g.,  uncertainties  associated  with  the  scanning  process  or  the  instrument).  POD 
determination  is  the  same  for  both  curves;  the  integration  of  the  signal  distribution  above 
the  SNR  threshold  (see  Eq.  4.1  and  Fig.  4.2).  Due  to  the  interaction  of  the  crack  and  the 
hole,  the  POD  values  will  not  be  equal  (notice  the  difference  in  the  amount  of  the 
distribution  above  the  threshold  for  each  case  in  Fig.  4.15)  but,  they  both  are  valid 
measures  of  POD  as  a  function  of  crack  length.  The  calculated  POD  curves  from  these 
two  cases  are  shown  in  Fig.  4.16.  The  difference  in  the  curves  results  from  the  hole/crack 
interaction,  not  a  difference  in  how  the  “signal”  was  defined,  and  shows  that  the  hole 
enhances  the  probability  of  detection.  Again,  the  crack  length  values  are  somewhat 
arbitrary,  since  the  distribution  was  chosen  for  illustrative  purposes  only  but  they  offer  a 
sense  of  scale  to  the  comparison. 
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Figure  4.16  Resulting  POD  curves  for  both  cases 


The  hole-with-crack  POD  curve  shows  that  a  hole  with  a  crack  has  a  smaller 
minimum  detectable  flaw  size  than  the  crack  only  for  a  given  probability  of  detection. 
This  is  directly  due  to  the  hole/crack  interaction.  The  effect  of  the  hole  increasing  the 
current  density  concentration  has  a  significant  effect  on  the  POD.  Even  though  the  POD 
curves  resulting  from  these  cases  can  be  compared  in  terms  of  detectability,  the 
hole/crack  interaction  does  not  allow  direct  comparison  by  use  of  a  baseline  hole  signal  to 
shift  one  curve  (vmless  you  completely  understand  this  interaction).  SQUID  performance 
will  have  to  be  determined  separately  for  each  kind  of  crack  type  (with  hole  or  without 
hole). 

Determination  of  POD  Curves 

Now  that  the  “signal”  has  been  defined  to  be  that  coming  from  the  hole  with 
crack,  POD  can  be  determined  from  the  hole-'with-crack  signal  distribution  (resulting 
from  the  imcertainties  in  the  system  parameters).  Sensitivity  analyses  will  be  done  on 
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selected  system  parameters  to  determine  how  they  affect  the  overall  signal  distribution. 
Monte  Carlo  simulations  using  this  signal  distribution  and  experimental  noise 
distributions  will  generate  realistic  POD  curves. 

System  Parameter  Uncertainty 

System  parameter  uncertainty  can  lead  to  intrinsic  noise  associated  with  the 
measurement  process.  The  uncertainty  in  a  parameter  manifests  itself  in  the  signal 
distribution.  To  determine  POD,  the  effect  of  the  uncertainty  in  each  system  parameter 
must  be  examined  to  determine  the  overall  signal  distribution.  The  BEM  measurement 
model  can  be  used  as  a  sensitivity  analysis  tool  to  look  at  how  the  magnetic  field  varies 
when  a  system  parameter  is  perturbed  within  a  range  associated  with  the  parameter’s 
uncertainty.  The  measurement  model  can  be  thought  of  as  representing  an  «-dimensional 
response  surface  where  n  represents  the  total  number  of  independent  parameters  used 
in  the  model  (e.g.,  liftoff,  crack  length,  and  pick-up  coil  diameter).  For  a  response  surface 
comprised  of  n  variables,  Xj,  such  that 

(4.12) 

7=1 

the  first  order  physical  sensitivity  factor,  PS^,  is  defined  as 

<5’'P  A'P 

PSj  =  (4.13) 

^Xj  Axj 

and  the  mean,  and  variance,  of  iPare  related  to  the  mean,  and  variance,  of 
each  of  the  variables  through  the  relations 
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(4.14) 


7=1 

By  looking  at  the  combined  effect  of  the  physical  sensitivities  (the  PS^s)  and  the 
distributional  characteristics  (the  a j  ’s)  of  the  system  variables,  along  with  some  insight 
obtained  by  understanding  the  underlying  physical  principles,  we  can  identify  factors  which 
influence  the  POD  distribution  most.  Those  system  parameters  that  have  a  relatively  large 
P5j  (sensitivity  analysis  will  show  a  steep  slope  on  the  response  surface)  or  cr,  (system 

parameter  distribution  will  be  spread  out)  will  be  most  important  in  the  determination  of 
POD  distributions  and  need  to  be  identified  either  through  analytical  or  experimental 
means.  It  may  be  necessary  to  treat  certain  aspects  of  the  inspection  process  as  “black 
boxes”  (e.g.,  data  acquisition  hardware)  and  examine  outputs  as  a  function  of  inputs 
without  mathematically  modeling  them.  Equation  (4.14)  requires  that  the  distribution 
characteristics  of  each  parameter  be  determined  or  estimated  in  order  to  define  the  overall 
signal  distribution  used  in  calculating  POD.  The  follovring  sections  will  examine  the 
parameters  of  scan  resolution,  plate  thickness,  pick-up  coil  liftoff  and  pick-up  coil  tilt 
angle.  The  uncertainties  associated  with  these  parameters  will  be  stated  in  terms  of  the 
error  (denoted  ABpp)  in  the  total  peak-to-peak  value  of  the  magnetic  field  5pp.  These 
distributions  are  all  determined  for  a  liftoff  at  or  near  3  mm  (minimum  distance  for  our 
SQUID  systems)  to  represent  the  maximum  capability.  For  other  liftoffs,  the  sensitivity 
analysis  would  have  to  be  done  again  to  generate  a  new  set  of  error  distributions. 
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Scan  resolution  uncertainty 

When  a  sample  is  scanned  under  the  SQUID,  the  user-defined  x-y  grid  of  the  scan 
determines  at  which  locations  data  will  be  acquired.  The  spatial  resolution  of  the  image 
is  very  dependent  on  the  fineness  of  the  grid.  It  is  common  practice  at  Vanderbilt  to  use  a 
1-mm  square  pixel  for  scanning.  Usually,  to  reduce  the  effects  of  noise,  multiple  data 
points  are  taken  within  the  l-mm  pixel  (in  the  scan  direction  only)  and  averaged  to  one 
value  for  that  pixel.  This  digitization  of  the  magnetic  map  can  lead  to  errors  in 
determining  the  signal  since  the  actual  dipole  extrema  are  not  likely  to  fall  exactly  on  a 
grid  point  (see  Fig.  4.17).  The  error  is  the  greatest  when  the  peak  is  in  the  middle  of  a 
pixel;  half  way  between  two  of  the  scan  lines  in  both  of  the  x  and  y  directions.  If  we  look 
at  the  variation  of  the  extrema  values  within  a  pixel  region,  the  sensitivity  of  the  signal, 
5pp,  to  scan  resolution  can  be  estimated. 


Figure  4.18  shows  the  maximum  error,  stated  as  a  percentage  of  the  total  peak-to- 
peak  value,  for  the  dipole  extrema  associated  with  a  5-mm  fastener  hole  with  a  4-mm 
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crack  at  a  liftoff  of  3  mm.  This  error  will  have  a  strong  dependence  on  varying  liftoff 
since  the  shape  of  the  dipolar  peaks  will  change.  The  error  will  also  be  dependent  on 
increasing  crack  length  but  to  a  lesser  extent  (see  shape  of  peaks  in  Fig.  4.10).  The 
dependence  of  scan  resolution  error  on  crack  length  is  a  subject  of  future  work.  The  error 
in  the  peak-to-peak  signal  ABpp  is  shown  for  a  1  mm  (scan  pixel  size)  variation  in  the  x 
andy  directions  (conducted  independently)  for  both  extrema.  The  maximum  error  for  the 
peak-to-peak  signal  is  the  sum  of  the  maximum  errors  for  each  extrema,  for  this  case, 
A5ppS  1.6%.  For  further  liftoff  distances,  the  peaks  become  less  sharp  which  would 
reduce  this  error  but  the  peak-to-peak  amplitude  decreases  with  a  greater  effect,  causing 
the  error  to  increase.  For  this  same  case,  but  with  the  liftoff  increased  to  6  mm,  A5ppS 
6%. 


1  mm 


Figure  4.18  Error  due  to  1  x  1  mm  scan  resolution 
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Since  it  is  equally  likely  that  the  dipole  extrema  will  fall  anywhere  within  the  scan  pixel, 
the  distribution  of  A5pp  associated  with  the  scan  resolution  uncertainty  can  be 
approximated  with  a  imiform  distribution  (in  terms  of  B^^),  =U(-0.016  j5pp,  0)  for  the 

3 -mm  liftoff  case.  Scan  resolution  imcertainty  always  reduces  Bp^,  which  explains  the 
minus  sign  in  the  distribution. 

One  comment  to  make  here  is  that  spatial  resolution  is  also  a  function  of  pickup 
coil  diameter  since  the  magnetic  field  is  being  integrated  over  the  face  of  the  pick-up  coil. 
In  other  words,  the  magnetic  flux  is  measured  and  represents  an  “averaged”  magnetic 
field  at  that  measurement  point.  A  rule  of  thumb  is  to  keep  the  scan  pixel  size  smaller 
than  the  pick-up  coil  diameter.  Pick-up  coil  flux  integration  was  done  for  all  calculations. 


Thin  plate  approximation  and  plate  thickness  uncertainty 

The  approximations  made  in  the  measurement  model  formulation  can  lead  to 
additional  uncertainties  in  the  calculated  magnetic  field.  One  approximation  was  made  in 
the  integration  through  the  thickness  of  the  plate  (see  Eq.  2.17,18).  It  was  assumed  that 
for  a  thin  plate,  the  1/r  term  was  approximately  constant  over  the  thickness  (no  z- 
dependence)  and  therefore  dQ  hds  with  the  thickness  h  assumed  constant.  This 
reduced  the  surface  integral  (Eq.  2.17)  to  a  contour  integral  (Eq.  2.18) 


[5  (c)]  =  ifo£  r  VF(^  ^  (g)  =  f  ^^^Qy^ds{Q)  (4.15) 

^  4;r  ■>«  r(c,0  ^  4;r  r(c,Q) 

If  we  examine  the  effects  of  this  assumption  more  closely,  the  limits  of  the  thin  plate 

approximation  can  be  shown.  The  error  arises  from  the  dependence  on  r,  which  varies 
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over  the  thickness  h,  as  shown  with  and  in  Fig.  4.19.  Only  r„pp,,  was  considered 
in  calculating  from  the  contour  integral,  which  was  then  multiplied  by  h  to  represent 
the  “2-D”  surface  integration. 


Figure  4.20  Integration  approximation  over  plate  thickness 


By  comparing  1/r  at  the  upper  and  lower  contours,  we  can  determine  the  magnitude  of  the 
error  associated  with  this  approximation 

W = (k  -  + (vc  -  yf  +  k  -  f" 
r,o.er  =  (k  -  ^y + k  -  yy + k  -  (z  -  h)y 

rewrite  in  terms  of  liftoff,  d  =  z^-z 


(4.16) 
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If  the  liftoff  is  much  greater  than  the  plate  thickness  {d  »  h),  the  ratio  is  approximately 
equal  to  one  and  the  error  would  be  small;  this  is  the  thin  plate  approximation.  However, 
for  increasing  plate  thickness  for  a  given  liftoff,  the  error  increases  (e.g. ,  doubled  at  h  =  d) 
and  can  no  longer  be  considered  negligible.  At  this  point,  the  plate  thickness  should  be 
subdivided  into  multiple  layers  (see  previous  Fig.  3.2)  and  the  integration  done  separately 
at  each  successive  layer,  then  added  together  (superimposed)  to  determine  the  total 
magnetic  field. 

The  plate  thickness  approximation  can  lead  to  error  in  the  calculated  magnetic 
field  but  this  can  be  treated  more  like  a  systematic  error  which  does  not  vary  if  the  scan 
geometry  (e.g.,  liftoff)  remains  the  same.  From  the  standpoint  of  scan-to-scan 
uncertainty,  it  is  the  variation  of  or  uncertainty  in  the  plate  thickness  h  which  is  of 
interest.  For  a  constant  value  of  injected  current  (measured  in  amps),  a  variation  in  h 
will  cause  the  injected  current  density  (measured  in  amps/m^)  to  change  and  thus  affect 
the  magnetic  field  amplitude.  The  magnetic  field  is  directly  proportional  to  the  injected 
current  density  which  is  proportional  to  the  plate  thickness  for  a  constant  injected  current 
and  plate  width 

=  =  l  +  (4.18) 

A  1%  variation  in  plate  thickness  will  cause  a  1%  variation  in  The  variation  could  be 
a  result  of  material  processing,  sample  fabrication,  or  error  when  measuring  the  thickness 
of  the  plate.  Estimating  the  uncertainty  associated  with  a  1-mm  plate  to  be  +  0.05  mm 
(10%  variation),  ABpp  associated  with  plate  thickness  uncertainty  can  be  approximated 
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with  a  normal  distribution  of  N(0, 0.025  B^^).  The  standard  deviation  of  this  distribution 
was  determined  by  assuming  that  the  maximum  error  (5%)  is  associated  with  the  2-sigma 
point  (approximately  95%)  of  the  A5pp  distribution. 

Pick-up  coil  liftoff  uncertainty 

The  magnetic  field  is  strongly  dependent  on  the  pick-up  coil  to  sample  distance 
(liftoff).  The  normalized  peak-to-peak  values  as  a  function  of  increasing  liftoff  distance 
for  a  3  mm  diameter  pick-up  coil  scanning  over  a  5  mm  hole  with  three  different  crack 
lengths  (2  mm,  3  mm,  and  4  mm)  are  shown  in  Fig.  4.20.  The  inset  graph  shows  a  detail 
of  the  curves  for  the  liftoff  range  of  3  to  4  mm.  Within  the  range  of  1  mm,  the  curves  for 
the  three  crack  lengths  are  approximately  linear  with  a  slope  of  30%/mm.  Normally,  the 
liftoff  distance  can  be  determined  within  +0.5  mm  and  therefore,  the  greatest  error  due  to 
uncertainty  in  liftoff  would  be  ±15%.  By  assuming  that  A5pp  due  to  liftoff  uncertainty  is 
normally  distributed,  the  approximate  distribution  for  A5pp  is  N(0, 0.075  Again,  the 
maximum  error  (15%)  is  used  as  the  2-sigma  point. 


62 


Figure  4.20  Sensitivity  of  to  liftoff  for  various  crack  lengths 


Pick-up  coil  tilt  angle  uncertainty 

The  BEM  measurement  model  has  assumed  that  only  the  z-component  of 
magnetic  field  (jBJ  is  being  measured  by  the  pick-up  coil.  If  the  pick-up  coil  is  not 
parallel  to  the  surface  of  the  plate,  will  be  reduced  and  the  other  components  of  the 
magnetic  field,  B^  and  B^,  will  contribute  to  the  signal.  These  components  are  associated 


with  the  i  and  j  terms  of  Eq.  (2.14) 
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Both  terms  have  a  common  multiplier  dependent  on  the  pick-up  coil  distance  to  the  top 
and  bottom  plate  surfaces  being  integrated  over.  The  effect  of  the  multiplier  on  each  of 
these  components  is  similar  to  the  previous  discussion  on  the  plate  thickness 
approximation  where  the  1/r  values  differed  at  the  upper  and  lower  contours  at  the  edge 
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of  the  plate.  Therefore,  by  similar  derivation,  the  magnitude  of  the  multiplier  term  is 
found  to  be  dependent  on  the  liftoff  to  plate  thickness  ratio. 
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This  multiplier  is  reduced  even  further  when  considering  the  geometric  term  associated 
with  the  coil  tilt  angle  components  in  the  x  and  y  directions.  A,  and  ,  shown  here  for 
the  x-direction. 
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It  can  be  seen  that  for  the  thin  plate  approximation  {h/d  «  1),  this  multiplier  term  is 
approximately  equal  to  zero  and  there  is  basically  no  magnetic  field  contributions  due  to 
and  By  For  plates  where  the  thin  plate  approximation  is  not  valid,  the  geometrical 
term  associated  with  the  coil  tilt  is  the  dominant  reduction  term.  For  example,  a  2  mm 
thick  plate  at  a  liftoff  of  3  mm  yields  O.lSsin  A^ .  For  small  tilt  angles  (<  5°),  the  B^ 
component  would  then  be  reduced  by  almost  two  orders  of  magnitude. 

These  two  components  appear  to  be  negligible  and  therefore  only  the  effect  of 
pick-up  coil  tilt  on  B^  will  be  considered.  The  peak-to-peak  value  of  5^  is  reduced  by  a 
geometrical  term  associated  with  the  tilt  angle. 


B 


pp 


=  5-_cosA 

liU=X  pp 


(4.21) 
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Tilt  angle  variability  is  small  in  the  laboratory  environment  since  the  SQUID  is  stationary 
and  samples  are  carefully  leveled  but,  in  the  field  this  variability  may  be  much  larger. 

For  now,  we  will  look  at  tilt  angles  <  5°,  which  result  in  errors  less  than  0.5%.  Since  the 
tilt  angle  will  always  result  in  a  reduction  of  the  distribution  for  A5pp  is  a  truncated 
distribution  (only  negative  values  used)  and  is  dealt  with  mathematically  by  taking 
negative  values  from  N(0,  0.0025 

Noise  Measurements 

Noise  distributions  must  be  characterized  in  order  to  determine  POD  and  POFA. 
Several  noise  distributions  have  been  extracted  from  existing  experimental  data  for 
different  SQUID  systems  and  measurement  techniques.  The  measurement  model 
simulates  dc-current  injection  and  experiments  using  direct  dc-current  injection  have 
large  noise  distributions  since  they  use  the  entire  baindwidth  of  the  SQUID  (dc  to 
-lOkHz)  and  therefore,  include  noise  over  these  frequencies  as  well.  Noise  conditions 
associated  with  other  experimental  techniques  are  more  representative  of  what  SQUID’s 
will  be  operated  in.  Techniques  based  upon  eddy  current  inducers  [29, 24]  use  lock-in 
amplifiers  at  a  particular  frequency  and  greatly  reduce  the  noise.  Figure  4.21  shows  a 
noise  time  series  extracted  from  experimental  data  (dc-current  injection  in  a  shielded 
environment)  by  taking  a  single  scan  line  away  from  any  known  flaw  in  the  specimen 
being  examined. 
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Time 

Figure  4.21  Experimental  noise  time  series 

As  can  be  seen,  the  noise  is  approximately  100  pT  peak-to-peak.  A  x^-test  confirmed  that 
this  noise  distribution  is  normally  distributed  at  a  5%  significance  level.  With  the 
inducing  plate  technique,  use  of  the  lock-in  amplifier  reduces  this  to  about  5  pT  peak-to- 
peak,  a  noise  reduction  factor  of  20.  Since  we  are  also  interested  in  taking  measurements 
in  an  unshielded  environment,  noise  distributions  have  been  extracted  from  scan  data 
taken  using  the  same  HTS-SQUID  that  was  used  in  some  of  the  experimental 
comparisons.  Figure  4.22  shows  noise  time  series  for  both  the  dc-measurement  (upper) 
and  the  inducer  measurement  using  the  lock-in  (lower).  The  dc-measurement  data  shows 
approximately  400  pT  peak-to-peak  noise  while  the  lock-in  measurement  data  shows 
about  2  pT,  a  noise  reduction  factor  of  200!  POD  values  determined  using  noise 
distributions  similar  to  that  shown  in  the  lower  part  of  Fig.  4.22  will  represent  present 
SQUID  capability.  Previous  tests  have  shovra  that  SQUID  experimental  noise  can  be 
assumed  to  be  gaussian  or  “white”  and  therefore,  the  noise  distribution  used  for  dc- 
measurements  was  N(0,  58.48  pT)  and  for  lock-in  measurements  was  N(0,  0.475  pT). 
These  noise  distributions  were  used  in  the  following  Monte  Carlo  simulation. 
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Noise 

(picoTesla) 


Figure  4.22  HTS-SQUID  noise  for  dc  (upper)  and  lock-in  (lower)  measurements 


Monte  Carlo  Simulation 


A  realistic  test  scenario 


Signal  strength  is  dependent  on  the  amount  of  dc-current  injected  into  the  sample. 


When  using  the  standard  dc-current  injection  technique,  it  is  a  common  approach  to 


increase  the  current  rmtil  adequate  SNR’s  are  achieved,  making  POD  primarily  dependent 
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on  how  large  a  current  source  one  might  have.  However,  the  dc  approach  is  not  practical 
and  is  not  likely  to  be  used  in  field  applications.  Instead  of  simulating  dc-current 
injection  directly,  the  BEM  measurement  model  will  be  used  to  approximate  the  fields 
associated  with  the  eddy-current  inducing  plate  [29]  mentioned  earlier  which  produces 
“dc-like”  currents  within  the  sample.  To  do  this,  an  estimate  of  the  current  density  must 
be  made  from  [30]  to  calculate  magnetic  fields  of  the  correct  magnitude.  The  inducing 
plate  normally  carries  a  current  of  approximately  0.04  mA/mm.  This  produces  eddy 
currents  on  the  order  of  0.01  mA/mm^  in  the  top  1  mm  of  an  aluminum  plate. 
Unfortunately,  for  the  thinner  plates  commonly  found  in  aircraft  fuselage  lap  joints,  the 
eddy  currents  in  the  bottom  of  the  plate  (in  the  opposite  direction)  reduce  the  overall 
signal  by  about  40%  for  a  1  mm  plate.  This  yields  a  final  estimate  of  current  density  of 
0.006  mA/mm^.  This  is  a  representative  value  associated  with  state-of-the-art  detection 
techniques  used  in  SQUID  NDE. 

The  test  sample  geometry  used  for  these  simulations  was  a  1-m^,  I -mm  thick 
aluminum  plate  containing  a  5-mm  fastener  hole  with  crack  lengths  ranging  up  to  4  mm. 
A  uniform  dc-current  of  0.006  mA/mm^  was  injected  in  a  direction  transverse  to  the  crack 
for  maximum  signal.  The  modeled  SQUID  system  utilized  a  3-cm  baseline  axial 
gradiometer  with  3-mm  diameter  pick-up  and  balance  coils.  The  plate  was  scanned  using 
1-mm^  x-y  pixels  at  a  pick-up  coil  to  plate  surface  liftoff  distance  of  3  mm. 

Steps  of  Monte  Carlo  distribution  sampling 

Monte  Carlo  simulation,  utilizing  the  results  from  the  BEM  measurement  model, 
was  used  to  sample  from  the  uncertainty  distributions  of  scan  resolution,  plate  thickness. 
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pick-up  coil  liftoff  and  pick-up  coil  tilt  angle  to  generate  an  overall  signal  distribution. 
The  resulting  signal  distribution  was  then  compared  with  the  noise  distributions 
associated  with  ac  and  dc  measurements  to  determine  POD.  A  95%  lower  confidence 
bound  was  then  generated  by  iterative  generation  of  multiple  signal  distributions  at  each 
crack  length.  The  Monte  Carlo  simulation  steps  are  outlined  here: 

1 .  For  each  crack  length,  sample  1000  times  from  each  of  the  four  system  parameter 
distributions  to  generate  1000  values  of  total  error  in  ->  (E  A5pp).  This,  added  to 
5pp,  constitutes  the  signal  distribution  to  be  compared  to  the  noise  distribution. 

-Distributions  used  related  to  system  parameter  uncertainty 
A5pp  (scan  resolution):  U(-0.016  0) 

A5pp  (plate  thickness):  N(0,  0.025  5pp) 

ABpp  (pick-up  coil  liftoff):  N(0,  0.075  B^^ 

A5pp  (pick-up  coil  tilt  angle):  negative  values  from  N(0,  0.0025  B^. 

2.  Within  the  same  sampling  loop  as  step  1,  sample  1000  times  from  the  noise 
distribution.  In  this  way,  each  B^  will  have  a  corresponding  noise  value. 

-Noise  Distributions 
dc:  N(0, 58.48  pT) 
ac(lock-in):  N(0,  0.475  pT) 

3.  Compare  each  signal-noise  pair  (1000  pairs)  to  determine  whether  B^^>  (noise  x 
SNR).  If  yes,  then  increment  “detected”  counter.  Otherwise,  increment  “missed” 
counter.  A  signal-to-noise  ratio  (SNR=2)  was  introduced  at  this  point  to  represent  a 
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detection  criteria  such  that,  when  the  signal-noise  pair  is  compared,  the  signal  value 
must  be  at  least  twice  as  large  as  the  noise  value  to  count  as  a  “detection”.  This  is 
basically  the  “threshold”  of  the  POD  integration. 

4.  Then,  for  that  particular  crack  length,  the  POD  =  (#  detected)/!  000. 

5.  For  determination  of  95%  lower  confidence  bound,  run  steps  1  through  4  100  times  to 
generate  a  distribution  of  POD  values.  Sort  the  100  values  in  ascending  order  (lowest 
POD  to  highest  POD).  The  95%  lower  confidence  bound  value  is  such  that  95%  of 
all  POD  values  in  this  distribution  are  above  this  bound;  this  corresponds  to  the  5"' 
POD  value  in  the  sorted  sequence. 

6.  For  each  crack  length,  run  steps  1  through  5.  Forty  crack  length  increments  were 
used  for  these  calculations.  For  the  dc-noise  calculation,  this  corresponded  to  a  range 
of  (0  -  4  mm).  For  the  ac-noise  calculation,  the  range  was  (0  -  0.05  mm). 

Using  the  sampling  numbers  listed  in  steps  1  through  6  resulted  in  a  total  number 
of  Monte  Carlo  simulations  of  (1000  /parameter)(5  parameters)(100  POD  loops)(40  crack 
length  increments)(2  noise  distributions)  =  40  million  points.  This  would  not  be  easily 
accomplished  experimentally;  at  one  measurement  per  minute,  it  would  take  over  75 
years! 

Monte  Carlo  generated  POD  curves 

Figure  4.23  shows  the  histogram  representation  of  the  signal  and  noise  data 
generated  (steps  1  and  2)  for  the  5-mm  hole  with  a  1-mm  crack  in  the  dc-measurement 
case.  The  noise  distribution  was  defined  to  be  one-sided  (a  rectified  normal  distribution) 
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since  the  signal  distribution  has  only  positive  values  (peak-to-peak  values  are  always 
positive).  In  this  way,  the  noise  distribution  is  represented  in  a  peak-to-peak  sense  so  that 
the  distributions  can  be  compared.  In  the  noise  histogram  (lower),  the  signal  histogram 
(upper)  has  been  appropriately  scaled  and  overlaid  (dashed  line)  to  show  the  region  of 
overlap  of  the  two  distributions.  It  is  this  region  which  is  of  interest  when  determining 
POD. 


Figure  4.23  Monte  Carlo  generated  signal  and  noise  distributions  for  the 
dc  case  (1  mm  crack) 
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Figure  4.24  is  the  dc-measurement  POD  curve  resulting  from  generating  similar 
distributions  at  each  crack  length.  The  abruptness  of  the  POD  curve  near  the  origin  is  due 
to  a  combination  of  the  relative  sharpness  of  the  signal  distribution  with  respect  to  the 
noise  distribution  and  the  one-sidedness  of  the  noise  distribution.  Also  shown 


Crack  Length  (mm) 

Figure  4.24  Ac/dc-measurement  POD  curves  showing  90/95  crack  lengths 
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in  this  figure  is  the  95%  confidence  bound  that  was  determined  from  steps  5  and  6.  This 
lower  bound  can  also  be  calculated  using  a  standard  one-sided  confidence  limit  formula 
[31].  For  this  case, 

,  (4.22) 

where  XpQi^  is  the  sample  mean  and  SpQ^  is  the  sample  standard  deviation  of  the  100 
Monte  Carlo  generated  POD  values.  Lower  confidence  bounds  using  Eq.  (4.22)  were 
much  narrower  (lower  confidence  bound  very  close  to  POD  curve)  than  those  determined 
through  direct  sampling.  It  was  decided  to  retain  the  Monte  Carlo  derived  confidence 
bounds  to  be  consistent  with  the  approach  and  this  would  also  give  the  “worst-case” 
lower  confidence  bound.  However,  Eq.  (4.22)  reveals  how  the  confidence  bound  depends 
on  the  number  of  data  points  or,  in  the  case  of  Monte  Carlo,  the  number  of  simulations 
and  therefore,  will  approach  the  POD  curve  itself  for  large  numbers  of  simulations.  The 
95%  lower  confidence  limit  applies  more  to  experimentally  based  POD  curves  that  are 
usually  generated  from  a  relatively  small  number  of  data  points.  For  simulated  data,  the 
lower  confidence  bound  can  be  made  to  basically  lie  on  the  POD  curve  if  enough  runs  are 
made  and  is  not  as  useful  as  a  concept  as  for  the  experimentally  derived  POD  curve. 

The  minimum  detectable  crack  lengths  corresponding  to  90%  probability  of 
detection  at  95%  confidence  are  1 .4  mm  for  the  dc-measurement  and  0.0134  mm  for  the 
ac-measurement.  The  very  small  minimum  detectable  crack  length  determined  for  the  ac- 
measurement  is  due  to  the  large  noise  reduction  through  use  of  the  lock-in  amplifier  in 
this  type  of  measurement. 
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POD  and  POFA 


In  the  Monte  Carlo  simulations,  the  signal  distribution  was  compared  directly 
with  the  noise  distribution  to  determine  POD.  For  direct  sampling  from  distributions,  the 
signal-to-noise  ratio  set  the  threshold  (a,^  in  standard  calculations  of  POD  (see  Eq.  (4.1)), 
which  set  the  detection  criteria.  To  examine  the  tradeoff  between  the  probability  of 
detection  and  the  probability  of  false  alarm  when  setting  a,;, ,  the  signal  and  noise 
distribution  data  generated  by  the  Monte  Carlo  simulations  can  be  presented  in  a  different 
format.  Refer  back  to  Fig.  4.3  and  Eqs.  (4.1)  and  (4.8)  regarding  the  description  of  POD 
and  POFA.  The  SNR  requirements  determine  where  d,^  is  to  be  placed,  thus  affecting 
POD  (larger  SNR  requirements  correspond  to  larger  minimum  detectable  crack  lengths). 

Figure  4.25  displays  the  dc-measurement  POFA  and  POD  curves  for  three  crack 
lengths,  including  the  minimum  detectable  crack  length,  as  a  function  of  threshold.  By 
setting  a  threshold,  a  value  for  POFA  and  POD  for  all  crack  lengths  is  determined.  As 
can  be  seen,  for  a  crack  length  of  0.5  mm,  the  90%  probability  of  detection  point 
corresponds  to  a  28%  POFA  (i.e.,  there  is  a  28%  chance  that  noise  will  be  mistaken  for  a 
signal).  There  is  only  a  single  POFA  curve  since  probability  of  false  alarm  is  determined 
from  the  noise  distribution,  which  for  these  simulations  is  constant.  Figure  4.26  shows 
the  ac-measurement  POFA  and  POD  as  a  function  of  threshold  for  the  previously 
determined  minimum  detectable  crack  length  of  0.01 34  mm.  In  both  the  dc  and  ac 
measurement  cases,  the  POFA  corresponding  to  the  minimum  detectable  crack  length  is 
approximately  zero.  This  is  because  the  determination  of  these  minimum  values  was  a 
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result  of  a  SNR=2  requirement  (equivalent  to  doubling  the  noise  distribution  in  Figs.  4.25 
and  4.26). 


Figure  4.25  Dc-measurement  POFA  and  POD  as  a  function  of  threshold 
for  three  values  of  crack  length 


Figure  4.26  Ac-measurement  POFA  and  POD  as  a  function  of  threshold 
for  the  minimum  detectable  crack  length 
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CHAPTER  V 


SUMMARY  AND  DISCUSSION 


Summary 

This  research  was  directed  towards  the  use  of  a  Superconducting  Quantum 
Interference  Device  (SQUID)  magnetometer  as  a  tool  for  nondestructive  evaluation  (NDE) 
to  detect  and  characterize  structural  damage  occurring  in  aging  aircraft.  The  primary 
advantage  of  using  SQUID’ s  in  NDE  over  other  techniques  is  their  ability  to  detect 
second  layer  damage  commonly  found  in  aircraft  structures.  A  probability  of  detection 
(POD)  analysis  is  being  done  to  validate  the  capability  of  SQUIDs  in  this  role.  A  POD 
analysis  can  be  done  using  real  experimental  data  or  simulated  experimental  data 
generated  from  a  measurement  model.  A  goal  in  this  research  has  been  to  develop  the 
latter  method,  where  simulation  represents  the  experimental  approach,  is  cheaper  and 
faster,  and  identifies  sources  of  unreliability  in  SQUID  NDE. 

The  approach  has  been  to  develop  a  measurement  model  simulating  the  scanning 
of  a  SQUID  over  a  dc-current  injected  sample  containing  a  crack.  The  modeled  crack 
and  fastener  hole  simulates  fatigue-cracking  conditions  found  in  aircraft  fuselage  lap 
joints.  Boimdary  integral  equations,  using  a  special  Green’s  function  incorporating  the 
crack,  are  used  to  solve  the  potential  problem.  This  special  formulation  eliminates  the 
need  to  model  the  crack  as  part  of  the  boundary  mesh.  From  the  BIE  potential  solution, 
the  magnetic  field  above  the  sample  can  be  calculated  through  numerical  integration  of 
the  Biot-Savart  law.  Special  techniques,  such  as  the  use  of  quadratic  interpolation 
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functions  and  a  quarter-point  formulation,  were  required  to  accommodate  the  singular 
behavior  at  the  crack  tip. 

Thin  copper  plate  samples  containing  a  thin  slot  simulating  a  crack  were  used  for 
experimental  comparisons  with  the  BEM  measurement  model.  Several  different  SQUID 
systems,  including  a  shielded  axial-gradiometer  SQUID  (MicroSQUID),  an  unshielded 
axial-gradiometer  SQUID  (SBIR  SQUID),  and  an  unshielded  planar-gradiometer  SQUID 
(Mobile  HTS  SQUID),  were  used.  Comparisons  of  scan  centerline  profiles  from  the 
experimental  magnetic  maps  show  good  agreement  (within  a  few  percent  in  the  region  of 
the  flaw)  with  the  measurement  model  for  all  SQUID  systems  and  test  conditions  used. 
Differences  can  be  attributed  to  sample-related  noise  (e.g.,  non-uniform  current  injection 
density). 

Monte  Carlo  simulation  utilizing  the  results  of  the  BEM  measurement  model  was 
used  to  generate  POD  curves.  The  uncertainty  distributions  resulting  from  a  sensitivity 
analyses  on  several  system  parameters  (scan  resolution,  plate  thickness,  pick-up  coil 
liftoff,  eind  pick-up  coil  tilt  angle)  were  combined  with  experimental  noise  distributions  to 
determine  POD  as  a  function  of  crack  length.  Threshold  analysis  resulted  in  comparative 
Probability  of  False  Alarm  (POFA)  and  POD  curves.  For  dc-measurements,  the 
minimum  detectable  crack  length  that  could  be  found  with  90%  probability  of  detection 
with  95%  confidence  was  1 .4  mm.  The  ac-measurements,  which  had  much  lower  noise, 
had  a  corresponding  value  of  0.0134  mm.  The  POFA  for  these  two  cases  was  negligible 
due  to  the  signal-to-noise  ratio  being  set  to  two  for  the  calculations.  These  small  crack 
lengths  suggest  that  additional  experimental  noise  factors  will  have  be  incorporated  into 
the  POD  analysis  before  realistic  SQUID  NDE  capability  can  be  accurately  quantified. 
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Discussion 


Factors  That  May  Decrease  POD  Capability 

The  relatively  small  values  of  minimum  detectable  crack  length  determined  from 
this  analysis  prompts  discussion  on  whether  this  reflects  the  true  capability  of  SQUIDs  in 
the  detection  of  real  cracks  near  fasteners  in  aircraft  lap  joints.  Experimental 
measurements  using  SQUIDs  on  fatigue  cracks  in  lap  joints  [36]  have  demonstrated  the 
difficulty  in  flnding  cracks,  adjacent  to  fasteners,  on  the  order  several  millimeters  in 
length.  For  the  POD  analysis  done  here,  we  have  considered  only  four  system  parameters 
(scan  resolution,  plate  thickness,  pickup  coil  liftoff  and  tilt  angle)  and  their  associated 
uncertainties.  Since  it  is  these  parameter  uncertainties  that  determine  the  overall  signal 
distribution,  and  therefore  POD,  it  is  apparent  that  significant  parameters  have  been  left 
out  of  this  analysis.  More  than  likely,  these  stem  from  sample-related  noise  associated 
with  real  measurements,  which  have  not  been  represented  in  either  the  measurement 
model  itself  or  through  the  parameter  uncertainty  distributions.  Sample-related  noise 
results  from  the  geometry  of  the  sample  or,  more  specifically,  how  that  geometry  affects 
the  current  flow  used  to  probe  the  sample.  Figure  5.1  shows  how  some  of  these  current 
paths  can  go  vmder  and/or  through  the  crack  and  the  fastener,  further  complicating  the 
detection  and  characterization  of  the  crack.  Some  of  these  issues  are  discussed  in  the 
following  sections. 
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Figure  5.1  Additional  current  paths  complicating  crack  detection 


Crack  shape 

Probability  of  detection  is  traditionally  stated  as  a  function  of  crack  length  to  be 
consistent  with  damage  tolerant  requirements,  which  are  stated  in  terms  of  crack  length. 
However,  crack  shape  is  an  important  issue  to  discuss  briefly  since  cracks  of  the  same 
“length”  but  different  shape  could  result  in  different  signals.  Currently,  the  BEM 
measurement  model  simulates  a  straight-through  crack  (the  same  length  through  the 
entire  thickness  of  the  plate)  but  it  is  possible  to  use  the  model  to  simulate  “3-D”  cracks 
to  examine  the  effect  of  crack  “shape”  on  magnetic  field.  By  creating  slices  of  varying 
shapes  and  depths  (see  Fig.  5.2),  a  pseudo-3D  crack  can  be  modeled.  The  total  signal  is 
then  the  superposition  or  sum  of  the  signals  calculated  at  each  of  these  slices.  However, 
as  discussed  in  Chapter  II,  the  model  does  not  allow  current  flow  perpendicular  to  the 
slices,  only  in  the  plane  of  each  layer,  and  cannot  model  current  going  under  cracks  and 
fastener  heads. 

For  this  work,  the  “signal”  used  in  the  POD  analysis  has  been  defined  as  the 


Figure  5.2  Pseudo  3-D  modeling  of  a  crack  using  slices 


peak-to-peak  value  of  the  magnetic  dipolar  signal.  The  dipolar  signal  has  characteristics 
that  are  related  to  physical  characteristics  of  the  crack.  The  amplitude  of  the  dipole  is 
proportional  to  crack  size.  More  specifically,  for  dc-current  injection,  the  signal  is 
proportional  to  the  cross  sectional  area  of  the  crack.  For  the  2-D  modeling  used  in  the 
BEM  measurement  model,  this  is  equivalently  crack  length.  It  is  also  known  that  the 
peak-to-peak  separation  is  proportional  to  the  liftoff  and  the  spreading  out  of  the  dipole 
signal  is  related  to  the  spatial  distribution  of  the  crack.  Since  crack  shape  affects  both  the 
amplitude  and  the  shape  of  the  magnetic  field,  then  this  pseudo  3-D  approach  would  be  a 
better  approximation  to  a  real  crack  that  has  shape  that  varies  with  depth. 

Ideal  cracks  versus  real  cracks  and  slots 

The  BEM  measurement  model  assumes  an  ideal,  perfectly  insulating,  infinitely 
thin  crack.  This  assumption  will  lead  to  differences  between  the  model  and 
measurements  made  on  samples  containing  either  wide  cracks  (slots)  or  real  fatigue 
cracks.  Most  fabricated  test  samples  used  in  our  laboratory  use  slots,  made  by  a  saw  or 
electrodischarge  machining  (EDM),  to  simulate  cracks.  A  combination  of  a  drilled  hole 
with  an  EDM  slot  is  an  approximation  to  a  crack  emanating  from  an  aircraft  fastener 


hole.  This  is  useful  in  preliminary  analyses,  especially  for  development  of  NDE 
techniques.  Fabrication  of  test  samples  made  this  way  is  simple  and  controllable, 
making  it  easy  to  build  a  test  set  representing  the  range  of  conditions  that  are  of  interest. 
But,  measurements  with  NDE  instruments  [18, 19]  have  shown  that  the  instrument 
response  from  a  slot  is  not  necessarily  the  same  as  that  from  a  fatigue  crack  of  the  same 
size  and  geometry.  Figure  5.3  schematically  shows  some  of  the  possible  variations,  each 
possibly  resulting  in  a  different  instrument  response. 


Figure  5.3  Slot  profiles:  (a)  wide  slot  (b)  narrow  slot  (c)  ideal  crack  (d)  real  crack 


More  realistic  crack  specimen  fabrication  can  be  accomplished  as  described  in  [20]: 

a)  Introduce  a  controlled  starter  notch  (usually  by  EDM)  into  the  specimen. 

b)  Promote  a  fatigue  crack  from  the  notch  by  cyclically  fatigue  loading  the  specimen 
with  axial  and/or  bending  loads. 

c)  Remove  starter  notch  to  retain  only  the  sharp  fatigue  crack  and  to  produce  a  desired 
specimen  surface  finish. 

However,  this  is  a  time  intensive  process  and  it  is  difficult  to  control  as  to  produce 
“standard”  cracks  of  various  lengths  used  in  system  characterization. 
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The  BEM  measurement  model  simulates  a  closed  crack  that  is  electrically 
insulated  along  the  entire  length.  We  have  shown  that  it  provides  a  reasonable 
approximation  to  a  thin  slot  for  liftoffs  appropriate  for  SQUIDs  commonly  used  in  the 
laboratory  (>  3  mm).  However,  for  real  fatigue  cracks,  it  is  possible  that  crack  closure 
may  cause  electrical  conductivity  across  parts  of  the  crack.  Probability  of  detection  will 
be  strongly  dependent  on  the  effects  of  crack  closure  since,  if  current  is  flowing  through 
the  crack  instead  of  around  it,  the  signal  (which  is  proportional  to  the  perturbation  of  the 
current)  will  be  greatly  reduced.  It  is  not  yet  understood  how  much  crack  closure  effects 
current  flow  since  it  is  likely  that  the  oxide  layers  that  form  on  the  crack  surfaces  will  act 
as  an  electrical  insulator  and  so,  only  those  contact  points  where  the  oxide  layer  is  not 
present  will  provide  conductive  paths. 

Additional  current  paths 

The  POD  analysis  has  been  based  on  simulated  measurements  of  a  crack 
emanating  from  a  fastener  hole  and  has  completely  ignored  the  effects  of  the  fastener 
itself.  This  is  most  likely  one  of  the  larger  sources  of  discrepancy  present  between  the 
model  and  real  measurement.  Since  most  fatigue  cracks  start  under  the  fastener  head, 
early  detection  is  difficult  and  usually  the  crack  has  to  propagate  beyond  the  fastener 
head  before  detection  occurs.  Another  critical  issue  is  that  of  contact  between  the 
fastener  and  the  hole  sides,  which  can  create  numerous  current  flow  paths  across  and 
through  the  fastener.  These  currents  may  be  too  difficult  to  model  directly  and  their 
associated  xmcertainty  might  have  to  be  represented  through  a  sample-related  noise 
distribution.  This  will  probably  have  to  be  determined  experimentally  by  examining 
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samples  with  holes  alone  and  with  fasteners  inserted  to  quantify  the  effect  on  the 
magnetic  signal. 

Geometry  factors 

Current  can  also  be  greatly  affected  by  the  geometry  of  the  surrounding  structure, 
particularly  edges  such  as  lap  joint  seams.  Figure  5.4  is  a  detail  from  Fig.  3.4  showing 
that  edges  can  produce  a  large  signal  amplitude  that  can  make  flaw  detection  difficult, 
especially  if  the  flaw  is  located  near  an  edge. 


Figure  5.4  Detail  from  Fig.  3.4  pointing  out  edge  and  flaw  signals 

Other  examples  of  edge  signals  can  be  seen  in  previously  shown  figures  (Fig.  1 .6  (b)  and 
Fig.  3.2).  Image  processing  techniques  can  partially  remove  the  background  slope  due  to 
these  edges  but  it  is  still  difficult  to  extract  the  signal  from  a  flaw  that  is  near  an  edge. 
Edges  and  other  geometrical  factors  (c.g.,  structural  support  members)  affect  the  ability 
of  the  SQUID  to  detect  a  crack  and  these  have  not  been  taken  into  accoimt  in  this  POD 
analysis. 
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Factors  That  May  Increase  POD  Capability 

The  peak-to-peak  value  of  the  magnetic  dipolar  map  is  the  present  definition  of 
“signal”  used  in  these  analyses.  Peak-to-peak  amplitude  is  a  very  limited  use  of  the 
information  available  in  the  mapping  since  it  reduces  the  entire  2-D  image  to  a  single 
value.  It  is  the  opinion  of  the  author  that  image  features  (e.g.,  shape  and  asymmetry)  of 
the  2-D  magnetic  field  map  are  just  as  important  as  the  peak-to-peak  signal  amplitude  in 
the  detection  and  characterization  of  cracks,  especially  near  fastener  holes.  For  example, 
if  we  look  at  the  magnetic  map  centerline  profiles  for  various  ideal  geometries  that  have 
the  same  “cross  sectional  length”  perpendicular  to  the  current  injection  direction,  we  can 
see  that  using  only  peak-to-peak  information  for  current  perpendicular  to  the  crack  can 
lead  to  detection  problems.  The  left-hand  side  of  Fig.  5.5  shows  the  profiles  for  the 
current  transverse  to  the  crack  at  three  different  liftoffs  for  a  9  mm  diameter  hole,  a  5  mm 
hole  with  a  4  mm  crack,  and  a  9  mm  crack.  The  signal  values  for  each  of  the  traces  are 
plotted  on  the  same  scale  and  represent  a  SQUID  system  using  3  mm  pickup/balance 
coils  with  a  3  cm  baseline.  From  this  information,  all  the  profiles  look  the  same  in  that 
they  are  dipolar  with  a  mostly  symmetrical  shape.  Based  on  peak-to-peak  amplitude  in  a 
single  direction,  we  can  not  distinguish  between  a  crack  alone,  hole  alone,  or  a  hole-with- 
crack  combination  unless  a  baseline  signal  amplitude  can  be  established  as  was  done 
earlier  in  this  POD  analysis  using  the  signal  associated  with  a  5  mm  hole  alone. 

However,  this  approach  would  not  work  here  since  the  5  mm  hole  is  not  a  constant 
feature  in  all  three  geometries. 

The  rotating  current  schemes  of  [30]  take  advantage  of  the  differing  two- 
dimensional  structure  of  these  three  geometries.  By  rotating  the  current  direction,  the 
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signals  associated  with  a  crack  will  go  through  a  maximum  (when  the  current  is 
perpendicular  to  the  crack)  and  a  minimum  (when  the  current  is  parallel  to  the  crack).  A 
hole  alone  will  not  show  this  cyclic  behavior  (unless  it  is  out-of-round).  On  the  right- 
hand  side  of  Fig.  5.5  are  the  corresponding  profiles  when  the  current  is  parallel  to  the 
crack  showing  the  9  mm  crack  having  a  flatline  (zero  signal)  and  the  5  mm  hole  with  a  4 
mm  crack  having  a  signal  corresponding  to  a  5  mm  hole  alone  and  the  9  mm  hole  is 
unchanged.  In  this  way,  the  three  geometries  could  be  distinguished  by  a  measure  of 
their  peak-to-peak  difference  at  the  two  current  injection  directions,  with  the  crack  alone 
being  the  largest,  the  hole-with-crack  combination  the  next  largest,  and  the  hole  alone 
having  no  difference.  Real  fatigue  cracks  will  have  some  small  signal  for  current  injected 
parallel  to  the  crack  (versus  zero  signal  for  the  ideal  cracks  of  Fig.  5.5)  but  will  still  show 
the  same  maximum/minimum  behavior  as  a  function  of  current  injection  direction  (the 
signal  difference  will  just  be  smaller  than  the  ideal  crack  case). 

It  may  also  be  possible  to  distinguish  between  these  three  cases  using  a  single 
current  direction.  By  looking  at  the  contour  maps  for  the  three  cases,  it  can  be  seen  that  it 
may  be  possible  to  use  dipole  extrema  eccentricity  as  an  identification  procedure.  Figure 
5.6  shows  the  contour  maps  for  these  cases  with  the  bidirectional  arrows  in  the  upper 
figure  representing  the  measures  related  to  extrema  eccentricity,  defined  here  as 


The  hole-with-crack  contour  plot  has  different  eccentricities  for  the  left  and  right  extrema 
of  the  dipole.  For  the  hole  alone  and  the  crack  alone,  there  is  left-right  symmetry  but  the 
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Centerline  position  (mm) 


Figure  5.5  Profiles  for  three  9-mm  “cross  section”  geometries  at  varying  liftoffs 
(Note:  all  profile  traces  are  on  the  same  scale) 


Figure  5.6  Contour  plots  showing  extrema  eccentricity 

eccentricity  is  different  between  the  two,  .  This  is  due  to  the  hole  having 

more  spatial  (physical)  extent  than  the  crack  in  the  direction  causing  the  signal  to  be 
more  stretched  out  in  comparison.  Since  this  feature  is  strongly  a  function  of  liftoff,  it  is 
difficult  to  determine  a  quantitative  measure  without  further  analysis. 

SQUID  NDE  systems  will  have  to  use  asymmetry  measures  (which  includes  the 
rotating  current  direction  technique)  for  crack  detection,  especially  for  cracks  originating 
fi:om  fastener  holes.  If  the  asymmetry  technique  is  used  only  to  find  the  maximum  and 
minimum  of  the  signal,  then  POD  may  still  be  determined  primarily  from  peak-to-peak 
amplitudes.  However,  if  additional  measures  of  asymmetry  (e.g.,  eccentricity  or  other 
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shape  factors)  are  to  be  used,  then  POD  analyses  must  incorporate  these  measures  into 
the  standard  peak-to-peak  analysis  to  reflect  the  true  capability  of  the  system.  This  is  an 
important  topic  to  be  addressed  for  future  development  of  SQUID  POD  methodologies. 

Other  Comments  and  Suggestions 

The  BEM  measurement  model  has  already  been  used  in  other  experimental  work 
involving  validation  [33,  34]  and  calibration  techniques  [35].  The  model  is  just  now 
beginning  to  be  used  for  POD  work,  which  was  what  it  was  originally  designed  for.  The 
measurement  model  can  provide  a  fast  and  accurate  way  to  simulate  scans  over  samples 
containing  injected  current,  making  it  a  very  useful  sensitivity  analysis  tool.  However,  if 
the  model  is  to  be  used  to  calculate  realistic  probability  of  detection  values  for  SQUID 
systems,  more  work  has  to  be  done  to  identify  those  sources  of  noise  which  are  causing 
the  large  discrepancy  between  experimental  measurements  on  real  fatigue  cracks  and  the 
detection  capability  determined  by  this  analysis.  Some  of  the  possible  sources  of  this 
noise  have  been  discussed  but  no  work  has  been  done  as  of  yet  to  quantify  their  effects. 

Also,  the  effect  of  the  crack  tip,  which  was  a  significant  model  formulation  issue, 
seems  to  be  somewhat  washed  out  at  the  liftoffs  presently  used  (>  3  mm).  It  may  be  that 
the  usefulness  of  the  BEM  measurement  model  may  not  be  fully  realized  until  the  liftoffs 
are  reduced  enough  to  see  the  large  magnetic  fields  in  proximity  of  the  crack  tip.  At 
smaller  liftoffs  (<1  mm),  using  less  accurate  models  (e.g.,  superposition)  will  lead  to  very 
large  errors.  Significant  improvements  in  reducing  liftoff,  requiring  a  very  short  distance 
between  superconducting  and  room  temperatures,  are  more  likely  to  be  achieved  in  high- 
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temperature  SQUID  systems.  The  measurement  model  and  the  POD  approach  are  readily 

applicable  to  these  types  of  systems  as  well. 

Some  additional  items  to  consider  for  future  development  include: 

-  automate  the  entire  process:  currently,  the  determination  of  POD  requires  several 
individual  computational  steps  and  it  should  be  possible  to  incorporate  all  of  these 
into  one  overall  program. 

-  refine  the  Monte  Carlo  simulation:  optimize  the  process  since,  as  more  uncertainty 
factors  are  included  in  the  sampling,  computational  time  will  increase  tremendously. 

It  may  be  possible  to  determine  functional  forms  for  the  POD  instead  of  using  Monte 
Carlo  methods  but  this  will  depend  on  the  distributions  associated  with  the  added 
uncertainty  factors.  Also,  we  need  to  define  how  to  compare  POD  curves  generated 
from  simulation  to  those  generated  from  experimental  measurements  (e.g.,  relevance 
of  lower  confidence  bound). 

-  examine  applicability  to  eddy-current  inducer  technique:  since  this  is  the  proposed 
technique  for  NDE  in  the  field,  we  need  to  determine  how  well  the  eddy  currents 
produced  by  the  inducer  sheet  can  be  modeled  by  the  dc-current  injection  of  this 
model. 

-  incorporate  signal  measures  other  than  Bz  peak-to-peak:  most  of  the  magnetic  map 
information  is  not  currently  used  (e.g.,  dipole  shape  characteristics).  Probability  of 
detection  may  be  enhanced  by  additional  measures  derived  from  this  unused 
information.  This  also  leads  into  the  use  of  vector  magnetometers  to  measure  and 
By,  which  has  yet  to  be  addressed  in  POD  determination. 
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APPENDIX  A: 


GREEN’S  FUNCTION  FORMULATION 
OF  LAPLACE’S  EQUATION 
FOR  ELECTROMAGNETIC 
CRACK  DETECTION 


Green’s  Function  Formulation  of  Laplace’s  Equation 
for  Electromagnetic  Crack  Detection 

T.  A.  Cruse,  A.  P.  Ewing,  and  J.  R  Wikswo^ 


1  Introduction 

The  following  note  summarizes  the  development  of  a  special  Green’s  function  for  crack  problems  in 
potential  theory.  The  formulation  is  for  a  single,  straight  crack  contained  within  or  intersecting  a 
regular  boundary.  The  formulation  of  the  special  Green’s  function  was  first  derived  for  the  elastic 
fracture  mechanics  problem.  The  current  formulation  is  a  subset  of  the  elasticity  problem  and 
is  formulated  as  a  Hilbert  problem  using  complex  variables.  The  formulation  is  validated  using 
boundary  element  modeling  of  the  Mode  III,  pure  (antiplane)  shear  fracture  mechanics  problem. 

The  Green’s  function  formulation  is  applied  to  modeling  the  magnetic  field  for  the  steady-state 
current  flow  through  a  two-dimensional  plate  with  a  plane  crack.  The  BIE  formulation  is  used  to 
obtain  the  boundary  integrals  necessary  to  compute  the  normal  component  of  the  magnetic  field  for 
an  arbitrary  remote  sensing  location.  The  boundary  integrals  involve  the  tangential  current  flow  on 
all  boundaries  including  that  for  the  crack.  The  explicit  form  for  each  of  the  two  crack  tip  singularity 
terms  is  derived. 


2  Singular  Potential  and  Green’s  Second  Identity 

The  field  problem  is  that  of  Laplace’s  equation  subject  to  mixed  boundary  conditions  which  may 
locally  be  Dirichlet  or  Neumann  conditions.  We  will  formulate  a  boundary  integral  equation  (BIE) 
for  this  problem  using  a  special  fundamental  solution  that  corresponds  to  a  surface  upon  which  the 
normal  derivative  of  the  solution  is  zero.  All  other  surfaces  will  have  user-specified  values  of  the 
potential  or  its  normal  derivative. 

The  field  equations  are  for  the  unknown  potential  4){p)  where  p(x)  is  the  Cartesian  field  point. 
The  BIE  will  be  formulated  from  a  Green’s  identity  using  the  fundamental  solution  of  Laplace’s 
equation  and  is  denoted  ip{p,  q)  where  p,  q  are  the  field  point  and  the  singular  point,  respectively. 
When  these  points  are  taken  to  the  boundary  they  will  be  denoted  as  upper-case  points  P,  Q.  Details 
supporting  this  note  are  given  in  [1]. 

The  fundamental  solution  in  two  dimensions  is  given  by  real  part  of  the  complex  logarithm 
function.  Recall  that  a  real  part  of  a  complex  function  is  one-half  of  the  function  plus  its  complex 
conjugate.  We  can  then  write  the  logarithmic  potential  as  the  following  real  form. 

V'(p,  q)  =  [log(2  -  c)  -h  log(f  -  c)]/2  =  7^1og(2:  -  c)  (1) 

where  c  =  Xp  +  iyp,  z  =  Xg  +  iyg  and  i  =  This  fundamental  solution  satisfies  the  following 

singular,  inhomogeneous  form  of  Laplace’s  equation. 

V^V’(P>  q)  =  2'ir'D(p,  q)  (2) 

In  this  equation,  V  is  the  Dirac  delta  function,  defined  by  the  relations  q)dV{q)  =  l(p)  and 

D(p,  ^)  =  0  so  long  as  q. 

The  integral  equation  formulation  for  potential  theory  derives  from  Green’s  second  identity, 
written  for  the  fundamental  solution  and  the  unknown  potential  function. 

[<^(g)V2^(p,  q)  -  iPip,  q)V^(l>{q)]  dV{q)  =  ^  _  ^(p,  Q)|^)  dS{Q)  (3) 

Fort  Flowers  Professor  of  Mechanical  Engineering,  Graduate  Student,  and  A.  B.  Learned  Professor  of  Physics, 
respectively 
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We  now  take  (j>{q)  to  satisfy  Laplace’s  equation  =  0  and  replace  the  Laplacian  of  the  fundamental 
solution  by  the  Dirac  function,  as  above.  The  result  is  an  integral  identity  for  the  values  of  (f){p)  in 
terms  of  the  totality  of  Dirichlet  and  Neumann  data. 

(4) 

3  Boundary  Integral  Equation  (BIE) 

The  above  identity  cannot  be  used  for  a  solution  to  the  unknown  potential  without  complete  knowl¬ 
edge  of  the  Dirichlet  and  Neumann  boundary  conditions.  In  general,  one  or  the  other  of  these  two 
boundary  conditions  is  specified  on  portions  of  the  boundary.  These  mixed  boundary  conditions  can 
be  written  as  follows. 


QeSv  (5) 

Q  e  Sn  (6) 

The  total  boundary  S  is  given  by  the  union  of  the  Dirichlet  boundary  So  and  the  Neumann  boundary 
Sn‘ 

The  Laplacian  operator  is  elliptic.  This  property  endows  the  potential  (j){p)  with  considerable 
smoothness.  The  boundary  values  of  0(Q)  will  be  taken  to  be  continuous  in  what  follows.  While 
the  derivatives  of  the  potential  are  also  continuous  on  the  interior,  the  normal  derivative  can  be 
discontinuous  on  the  surface.  That  fact  is  not  central,  however,  to  the  following  developments. 

Under  these  circumstances,  the  potential  (l>{Q)  and  its  normal  derivatives  are  continuous  at  the 
boundary.  We  can  now  subtract  the  boundary  value  of  the  potential  at  an  arbitrary  boundary  point 
F  from  the  boundary  integral,  and  add  it  back  in  as  follows. 

The  free  term  multiplying  the  potential  at  P  can  be  integrated  in  closed  form  for  all  interior  points 
p  and  the  result  is  independent  of  the  actual  surface  shape,  so  long  as  the  surface  is  closed. 

(8) 

We  now  obtain  an  equivalent  Green’s  identity  for  which  all  integrals  are  continuous  as  p  — >  P. 

2a<hp)  -  -  fm)  -  (») 

Green’s  identity  can  now  be  taken  to  the  boundary,  letting  p  — ►  P,  in  a  smooth  manner.  The 
resulting  boundary  integral  equation  (BIE)  for  the  potential  theory  formulation  is  given  by  the 
following  relationship. 

0 = jm)  -  - 1  m  do) 

The  BIE  is  a  constraint  equation  that  holds  for  all  harmonic  functions.  The  defined  constraint  is 
between  the  Dirichlet  and  Neumann  boundary  conditions  for  the  problem.  All  well-posed  problems 
can  be  reduced  to  a  well-defined  integral  equation  which  can  be  solved  for  all  unknown  boundary 
conditions.  Following  the  solution  of  the  BIE,  one  can  substitute  the  totality  of  the  derived  Dirichlet 
and  Neumann  boundary  data  into  the  Green’s  identity  for  the  interior  potential  function.  This 
constitutes  the  full  solution  of  the  boundary  value  potential  theory  problem. 


m = f{Q) 

d<f>\ 


dn 


=  9{Q) 
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4  Green’s  Function  for  Crack  Problems 


We  now  seek  to  derive  the  special  fundamental  solution  for  two  dimensional  crack  problems.  For 
the  EM  problems  we  will  be  using  the  potential  function  to  represent  the  current  potential  for  the 
steady-state  problem.  The  derivative  of  the  potential  normal  to  the  crack  surface  will  be  zero,  for 
the  assumption  of  a  perfectly  insulating  crack  surface.  We  seek  a  Green’s  function  for  the  resulting 
Neumann  potential  problem  of  the  insulating  crack.  That  is,  we  seek  a  function  whose  derivative 
normal  to  the  plane  crack  is  zero  at  the  crack  surfaces.  In  the  steady-state  EM  problem,  the  crack 
is  taken  to  be  an  infinitely  thin  insulating  surface.  We  define  the  surface  to  be  F.  Further,  we  will 
take  r  to  be  a  single  straight  line  parallel  to  the  x-axis  in  the  two  dimensional  problem.  The  source 
of  the  derivation  we  will  use  is  [2],  [3]. 

The  new  fundamental  solution  or  Green’s  function  is  the  combination  of  two  functions.  The 
first  is  the  same  logarithmic  potential  we  started  with;  the  second  is  a  function  which  will  cancel 
the  values  of  the  y-derivative  of  the  logarithmic  potential  —  at  the  crack  surface.  That  is,  the 
y-derivative  of  the  sum  of  the  two  potentials  is  to  be  zero  at  the  crack  surface.  We  achieve  this  by 
finding  the  solution  to  the  boundary  value  problem  for  the  second  potential  problem  that  cancels 
the  values  of  the  y-derivative  of  the  logarithmic  potential,  at  the  crack. 

The  problem  whose  solution  we  now  seek  is  called  a  Hilbert  problem  [4,  5].  The  formulation 
consists  of  finding  the  potential  for  an  infinite  region  with  a  crack  surface  F  on  which  we  have 
specified  boundary  conditions.  The  boundary  conditions  for  the  Hilbert  problem  are  to  be  the 
Neumann  conditions  involving  the  y-derivative  of  the  log-potential  at  the  crack.  For  superposition, 
the  resulting  derivative  is  the  negative  of  the  boundary  values  for  the  Hilbert  problem  and  are  taken 
to  be  P{t)  for  t  €  F. 

We  now  derive  the  complex  potential  h{z)  whose  boundary  values  on  F  are  P{t).  The  boundary 
values  for  our  problem  are  taken  to  be  the  same  on  both  the  upper  and  lower  side  of  the  crack  F. 

h+{t)=p{t) 

h-{t)  =  Pit)  (11) 

The  following  combinations  of  these  boundary  conditions  hold  for  the  Hilbert  problem  we  have 
defined.  The  general  solution  allows  the  potential  to  be  discontinuous  at  the  surface  F. 


=  0 

h-^{t)  +  h‘-{t)  =  2/3{t) 


(12) 


The  solution  for  this  Hilbert  problem  is  given  as  follows  [2]. 


1  f  13 ,  P{z) 
'Kyfz^  —  0?  Jr  t  —  z  yj  —  o? 


(13) 


In  this  equation,  the  integral  is  performed  on  the  positive  side  of  the  crack,  from  the  left  hand 
end  t  =  —a  to  the  right  hand  end  t  =  H-a.  The  function  of  integration,  P{z)  is  a  positive  power 
polynomial  which  can  be  determined  by  the  behavior  of  h{z)  at  infinity.  These  conditions  require 
that  P{z)  =  0.  A  discussion  of  the  complex  variable  integration  method  for  this  integral  is  given  in 
[5].  In  what  follows,  we  will  compute  only  the  real  part  of  this  integral. 

We  can  now  take  the  boundary  values  /3(f)  to  be  such  that  they  cancel  the  boundary  conditions 
of  the  y-derivatives  of  the  logarithmic  fundamental  solution.  This  superposition  of  two  functions 
gives  us  a  new  fundamental  solution  whose  y-derivative  is  zero  on  the  crack  surfaces.  Recalling  the 
fundamental  solution  to  be  given  by  the  following  expression 


i)[z,c) 


[log(2r  -  c)] 

i  [log{z  -  c)  +  \og{z  -  c)] 


(14) 
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we  can  derive  the  y-derivative  of  the  logarithmic  potential  at  the  crack,  which  is  parallel  to  the  x- 
axis.  We  will  take  the  negative  of  the  derivative  as  our  boundary  conditions  for  the  Hilbert  problem. 
The  final  solution  we  seek  will  have  zero  y-derivative  at  the  crack  by  adding  the  Hilbert  problem 
potential  to  the  fundamental  solution  potential. 


d’ip{z,c) 

1 

i 

i 

dy 

z=x"2 

X  —  c 

X  —  c 

=  -0{x) 


(15) 


It  is  seen  that  the  boundary  condition  is  a  real  function  on  the  crack  surface  F. 

The  Green’s  function  is  to  be  a  real  function  with  the  boundary  conditions  satisfying  the  values 
13 (x);  we  will  take  this  real  function  to  be  given  by  G{p,  q;  a).  The  Hilbert  problem  potential  is  then 
given  by  the  following  integral. 


G{p,q\  a)  =  -^71 


1  I"  y/o?  —  <2  ^  i  i  ^ 

—  Ji-  t  —  z  \f-c  t-cj 


dt 


(16) 


The  above  integrals  for  the  Green’s  function  terms  exist  in  closed  form,  as  can  be  seen  in  [5].  We 
need  only  the  real  part  for  each  of  the  above  terms,  as  in  [3],  given  as  follows 

G{jp,  q;  a)  =  Tllhiz)] 


=  -2^ 

where  the  integral  I{z)  is  defined  below 


i _ //(2)  -  j(c)  ijz)  -  /(g) 


V 


z  —  c 


)] 


I{z) 


-I 

TT  J-a 


dt 


_  ^2  _  2-j 


(17) 


(18) 


for  all  ^  except  along  the  crack  where  I{z)  has  the  boundary  values 


7+(i)  = 

r{t)  = 


t 


—  0? 
—i\/t^  -  0?  —  tj 


ter 


(19) 


The  complex  potential  h{z^  c)  has  been  used  to  find  a  Hilbert  problem  solution  for  the  boundary 
data  derived  from  the  y-derivatives  of  the  fundamental  solution.  These  terms  will  be  used  for  the 
evaluation  of  the  BIE  terms  involving  the  y-derivative. 

The  BIE  also  requires  that  we  have  the  fundamental  solution  potential,  and  not  just  its  y- 
derivative.  Therefore,  we  have  to  integrate  the  function  h{z)  with  respect  to  zji  to  have  the  function 
whose  boundary  conditions  on  the  crack  we  have  just  used.  The  integration  result  is  given  as  follows 
[6] 

0(p>9;o)  ='7^  l^T  J  h{z,c)dz  =  -^n[J{z,C]a)  -  J{z,c-,a)] 
where  c)  is  given  by  the  following  result 


(20) 


such  that 


J(z,c;a)  =  log  2 

dz 


\/z^  —  o^y/c?  —  (^  +  cz  — 

(z  +  \/z2  —  a2)  (c  +  \/c^  —  a^) 


(21) 


(22) 


{z  —  c)y/z^  — 

We  now  have  the  two  Hilbert  problem  forms  that  are  needed  to  modify  the  BIE  for  the  crack 
problem.  The  new  potential  J(z,c;a)  depends  on  both  field  points  as  well  as  the  crack  size.  The 
crack  for  this  formulation  is  limited  to  a  straight  crack  along  the  real  axis  with  its  center  at  the 
origin.  For  computational  purposes,  we  use  a  simple  shift  and  rotation  of  the  coordinate  axes  to  put 
the  crack  at  any  user-defined  location  and  orientation  relative  to  the  physical  geometry. 
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5  BIE  with  the  Special  Green’s  Function 


The  final  formulation  for  the  crack  problem  is  obtained  by  taking  the  sum  of  the  original  fundamental 
solution  with  the  solution  of  the  Hilbert  problem.  This  modified  fundamental  solution  which  satisfies 
the  insulated  condition  on  the  crack  surface  is  denoted  a  Green’s  function  for  the  crack.  That  is, 
the  Green’s  function  is  an  harmonic  function  containing  the  fundamental  singularity  property  and  it 
satisfies  the  boundary  conditions  on  the  crack.  This  fundamental  solution  is  denoted  by  q)  and 
is  given  by  the  real  part  (denoted  by  the  TZ  notation)  of  the  two  terms.  The  original  fundamental 
solution  term  has  been  normalized  by  dividing  it  by  27r. 


q]  a) 


_1_ 

27r 

2 

47r 


7^1og(2  -  c)  -  c;a)  -  J{z,c\a)} 


TZ  [log(2:  —  c)  +  log(2r  -  c)  -  J{z,  c;  a)  +  J(z,  c;  a)] 


(23) 


The  modified  fundamental  solution,  or  Green’s  function,  for  the  crack  is  shown  containing  the  crack 
length  as  a  parameter.  This  is  to  reinforce  that  the  presence  of  the  crack  is  now  embedded  in  the 
fundamental  solution  terms.  The  Green’s  function  derivative  in  the  normal  direction  at  the  surface 
point  Q  will  also  be  needed  for  the  BIE  formulation.  That  result  is  now  given. 


d^{p,Q-,a) 

II 

Tijjj  “f”  iTly 

Tlx  ~~ 

dn 

Z  —  C 

H - 'j: . . 

z  —  c 

1  ^ 

'Ux  +  lUy 

z  —  c 

1  1  _  r 

27r 

1 

1 

-iny)  / J 

^  V 


z  —  c 


z  —  c 


{Ux  +  iUy)  ( I{z)  -  I{c)  I{z)  -  I{c) 


z  -  c 


z  —  c 


(24) 


-  o2 

The  BIE  identity  can  now  be  written  on  the  total  of  the  regular  surface  and  the  crack  surface. 

.WQio) 


0  = 


dn{Q) 

d<P{Q) 


dS{Q)  + 
dS{Q) 


(25) 


The  effect  of  the  special  Green’s  function  formulation  is  now  seen.  Both  integrals  on  the  crack 
surface  T  are  zero.  The  first  is  zero  due  to  the  Green’s  function  having  zero  normal  derivative  on  the 
crack,  and  the  second  due  to  the  insulating  boundary  conditions  for  the  unknown  potential  function. 

Thus,  we  taJce  as  the  final  BIE  for  the  insulated  crack  potential  theory  problem  to  be  the  following 
equation. 

0 = Jm)  -  <“) 

The  BIE  models  the  crack  through  the  Green’s  function,  and  not  through  the  boundary  of  the 
crack.  The  form  that  has  been  derived  places  the  crack  at  the  origin  and  oriented  along  the  x-axis. 
However,  the  existing  computer  program  allows  the  user  to  specify  the  global  location  of  the  crack 
and  its  orientation.  The  crack  can  be  contained  within  the  region  or  can  intersect  one  or  more 
physical  surfaces.  Again,  the  BIE  code  has  been  written  to  accept  these  cases. 

Thus,  the  BIE  for  the  insulated  crack  problem  contains  only  the  uncracked  surface  in  the  BIE. 
Application  of  the  special  Green’s  function  to  the  elasticity  problem  has  demonstrated  the  very 
high  accuracy  of  this  formulation  for  crack  problems,  as  in  [6].  Further,  applications  of  the  special 
Green’s  function  formulation  have  demonstrated  a  high  degree  of  accuracy  for  cracks  at  holes  and 
other  stress  concentrating  geometries  [7]. 
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6  Interior  Derivatives  of  the  Potential 


The  interior  potential  for  the  Green’s  function  formulation  is  given  by  the  identity 


dn 


dS{Q) 


(27) 


The  derivatives  in  the  x-direction  and  the  y-direction  of  the  interior  potential  are  needed  for  the 
full  electromagnetic  field  formulation. 


d4>{p) 

d{xc,  Vc) 


dn{Q)d{xc,yc) 


Is 


d'i{p,Q;a)  d<f> 
d{xc,yc)  dn 


Q 


dS{Q) 


(28) 


The  derivative  at  the  point  c  of  the  kernels  given  in  Eq.  23  and  Eq.  24  are  required  in  order  to 
compute  the  gradient  of  the  interior  potential.  These  kernel  derivatives  are  given  as  follows. 


9^(p,  Q]  a) 
d(xc,  yc) 


(1,-f)  fIiz)-I{c)V 
\/c2  —  V  z  -  c  ) 


(29) 


The  identity 


dc  (c  -  2)Vc2  -  a2  (z  -  c)Vc^  -  a2  ^  ^ 

has  been  applied  using  the  symmetry  of  J{z,  c;  a)  with  respect  to  z,  c. 

The  derivative  at  c  of  the  normal  derivative  kernel  function  in  Eq.  24  is  given  in  the  following 
set  of  terms. 


a^^(p,Q;a) 
dnd{xc,  yc) 


(ux  +my)(l,i)' 

{z  -  c)2  _ 

{nx  +  iny){l,i)  (  d  I{z)  -  I{c)\ 
y/z^  —  a?  \9c  z-c  ) 
(n^4-m„)(l,-i)  fdl(z)-l{c)\ 
s/z"^  —  0?  \  5c  2  —  c  / 


The  remaining  derivatives  with  respect  to  z  above  eire  given  in  the  following  result. 


5J(2)-J(c)  J(c)  J(2) 

dc  z-c  (2  -  c)V'c2  -  a2 


(31) 


(32) 


7  Crack  Tip  Field  Intensity  Factory 


One  of  the  key  parameters  in  crack  problems  is  the  so-called  crack  tip  intensity  factor  (GIF  =  K) 
derived  from  the  derivatives  of  the  potential  at  each  crack  tip.  The  GIF  is  given  by  the  following 
limit 

(33) 

for  z  ^  c.  The  intensity  factor  can  be  derived  for  any  direction  of  approach  to  either  of  the  crack 
tips.  For  convenience,  we  will  take  the  point  c  =  ±(a  H-  r)  where  the  real  variable  r  is  the  distance 
from  the  crack  tips  and  is  taken  to  be  along  the  x-axis.  Thus,  c  =  c  for  the  limiting  processes.  The 
GIF  is  then  given  by  the  following  integral  identity  from  Eq.  4. 


^{x,y)\±a 


=  lim  y/2'Kr 
r-»0 


d{x,y) 


=  lim 

r—*0 


d^^{p,Q) 
d{x,  y)dn{Q) 


a^(p,g)5.^(Q)y 

d{x,y)  dn{Q)J^^ 


(34) 
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Carrying  out  the  limits  on  the  two  kernel  functions,  one  finds  that  the  derivative  d(p/dx  is  zero 
ahead  of  the  crack.  That  is,  the  gradient  of  the  flux  parallel  to  the  crack  is  zero  for  potential  theory. 
The  singular  term  is  contained  solely  in  d(j)/dy^  the  gradient  transverse  to  the  crack,  ahead  of  the 
crack.  One  can  also  find  from  a  study  of  the  kernels  that  the  gradients  for  points  on  the  crack 
surfaces  are  swapped  and  the  sign  changed.  Thus,  for  the  top  crack  surface  d(t)/dx  is  singular  with 
the  GIF  equal  to  minus  that  of  the  approach  ahead  of  the  crack,  output  in  the  code  as  K,  The 
derivative  d(f)/dy  is,  by  the  Green’s  function,  zero  along  the  crack  surfaces. 

The  two  kernels  that  are  needed  in  order  to  compute  K  for  the  limit  using  points  ahead  of  the 
crack  are  given  as  follows,  using  the  right  hand  crack  tip  c  =  +a.  The  code  allows  the  negative  crack 
tip  computation  as  well. 


/_52^(p,Q) 

r-^0  ay 


47r 

1 


7^ 


{{ria:  -\-iny) 


y/Aira 


n 


{z  —  a)y/z^  —  o? 
i[Iiz)-I{a)\ 


z  —  a 


(35) 

(36) 


8  Boundary  Element  Algorithm 


8.1  Quadratic  Straight  Line  Segments 

The  boundary  element  method  (BEM)  replaces  the  actual  boundary  with  an  approximate  boundary 
discretized  as  a  finite  set  of  simple  elements.  The  current  algorithm  uses  straight  line  segments 
in  order  that  we  can  implement  analytical  integrations  of  the  kernels  in  Eq.  23  and  Eq.  24.  The 
boundary  functions  are  then  represented  by  quadratic  interpolations  over  each  boundary  segment. 
Continuity  of  the  potential  is  enforced,  but  not  continuity  of  the  flux,  between  the  boundary  elements. 
The  analytical  integrations  of  the  two  kernels  times  quadratic  polynomial  interpolations  of  the 
boundary  data  are  computed  on  an  element-by-element  basis. 

The  quadratic  boundary  data  model  replaces  each  of  the  boundary  functions  with  the  following 
general  form  of  inter polative. 

=  Ao  +  Aiz  +  A2z‘^  (37) 

The  current  BEM  algorithm  is  based  on  boundary  elements  defined  by  end  nodes  and  mid-length 
nodes  for  each  boundary  element,  A5i  =  Az/B,  The  general  form  in  Eq.  37  can  be  replaced  by  a 
form  specifically  tailored  to  the  element-by-element  integration  for  the  BEM  algorithm.  We  define 
the  endpoint  nodal  values  of  each  interpolative  by  Fi,  F2  and  the  midside  nodal  value  by  F3,  The 
general  quadratic  interpolative  is  then  given  by 


F{z)  = 


F\  —  F2 

Az 


-ZQ 


Fi 


Az 


■^2^  .  0-^1  +-^2^2 

z  +  ^  /  A  Zq 


{Az^ 


Fi  +  F2  ,  0-^1  ^2  2 

-4  ..  .A  ZQZ  +  2  , .  ,A  -g 


(A2)2 


{Azf 


+  F3 


1~4 


+  8 


-4- 


{Azy  '  (Az)2  ^(A;^) 


(38) 


The  midside  node  for  each  segment  is  denoted  as  zq. 


8.2  Integration  of  the  Boundary  Terms 

The  integrals  for  each  boundary  element,  with  the  complex  variables  used  herein,  are  converted 
to  complex  integrals  with  respect  to  dz  —  BdS  where  B  =  iux  —  ny.  It  can  be  noted  that  A  = 
Tlx  +  iuy  =  —iB  in  Eq.  24.  The  boundary  integrals  are  given  symbolically  in  the  following  forms. 

J(c)  =  ^  F{z)  c;  a),  J  (39) 
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Substitution  of  the  quadratic  interpolative  from  Eq.  38  into  the  above  integral  results  in  the  following 
forms  of  boundary  element  terms  for  the  end-point  nodes  {N  =  1, 2)  and  the  midside  node  (N  -3). 

I{c)  =  A/o-(-l)^(A/i  -  A/2)- A/3  +  A/4  for 

=  A/5  +  A/s  +  A/7  for  JV  =  3  (40) 

The  terms  /q  and  h  result  from  integration  of  the  kernel  functions  times  the  constant  interpolation 
terms  in  Eq.  38.  The  terms  I2  and  I3  are  for  the  linear  terms  in  the  interpolation,  and  the  term  I4 
is  for  the  quadratic  interpolation  of  the  kernels.  The  terms  /j  and  h  change  sign  for  the  two  end 
nodes  in  Eq.  38.  In  the  case  of  the  integrals  relative  to  the  midside  node,  the  term  /s  is  for  the 
constant  term  in  the  interpolative,  the  term  /e  is  for  the  linear  term,  and  /y  is  for  the  quadratic 
term. 

The  kernel  functions  for  the  integrations  are  repeated  in  a  convenient  form  below. 


4'(p,  g;o) 
q]  a) 

dn 


^n[2log{z  -  c) 
47r  ^  -  c 


-  J{z,c;a)  +  J{z,c]  a)] 

,  .  .dJ{z,c)  . 

-  (rii  +  tny) — — - h  (ux  +  triy) 


dJ(z,cy 

dz 


(41) 


Substituting  the  identity  dz  =  BdS  in  the  above  equation,  and  integrating  over  the  segment  A5i, 
we  have  the  following  integrals  to  evaluate.  Since  we  are  using  straight  line  segments  in  the  BEM 
algorithm,  the  term  B  is  a  constant  over  each  segment. 


)^(p,  q\  a){p,  q;  a)dS 


,2, 2 


2'g^(p,  q-,a) 

’  dn 


dS 


4^  Ia  [21og(2  -  c)  -  J{z,c;a)  +  J{z,c\a)\dz 


dJ{z,c)  dJ{z,c) 
dz  dz 


dz  (42) 


The  final  results  are  obtained  by  taking  only  the  real  part  of  the  above  integrals. 

We  will  first  consider  the  integrals  of  the  potentials  of  the  uncracked  body.  The  integrated 
kernel  functions  are  labeled  AUq,  ATq,  AU\,  etc.,  referring  to  the  constant,  linear,  and  quadratic 
interpolatives.  The  following  summarizes  those  integration  results. 


27rAt/o 

27rATo 

2nAUi 

27rATi 

27rAC/2 

27rA!/2 
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—t 
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(43) 


The  following  terms  are  the  integrals  of  the  special  Green’s  function  terms.  The  integrations  will  be 
given  in  terms  of  the  source  point  c  as  the  conjugate  results  are  obtained  by  simple  substitution  of 
c  for  c. 


47rBA17o  =  -\{z  -  c)J{z,c-,a)  -  z  +  +  I{c)log{z  ^  +  .  ■  {z,c)] 

At: ATq  =  i[J[z,  c;  a)  -  J{z,  c;  a)] 
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(44) 


Again,  we  recall  that  the  BEM  equations  are  obtained  by  taking  only  the  real  part  of  the  above 
integration  results.  A  linear  system  of  equations  is  generated  by  the  BEM  algorithm.  The  system 
of  equations  is  automatically  swapped  during  assembly  such  that  the  unknown  boundary  conditions 
appear  in  a  single  vector.  The  system  of  equations  is  square  so  long  as  there  is  only  one  unknown 
flux  value  at  each  boundary  node. 


8.3  Integration  of  the  Interior  Gradient  Terms 

The  previous  integrations  are  used  to  compute  the  boundary  integrals  for  the  BIE  as  well  as  for 
the  evaluation  of  the  interior  potential  function.  No  we  consider  the  evaluation  of  the  derivatives  of 
the  interior  potential  function.  We  will  take  the  gradient  of  the  potential  at  an  arbitrary  interior 
point.  The  integrations  of  the  kernels  in  Eq.  29  and  31  will  again  be  broken  up  into  the  terms  for 
the  uncracked  potential  and  the  terms  for  the  special  Green’s  function.  The  integration  results  for 
the  interpolated  boundary  functions  times  the  kernels  will  again  be  denoted  by  Uq^  Tq,  Ui,  etc. 

The  uncracked  potential  function  integrals  are  given  as  follows. 


^t^(0,l,2) 


(9^(p,  Q;a)  dz 
d{xc,  Vc)  B 


(45) 
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uncracked  potential  can  be  integrated  from  each  source 


2nUo 

2tvUi 

2ttU2 

The  above  three  results  use  the  fact  that  the 
point  as  c  =  0  without  loss  of  generality. 

27rTo  = 

27rTi  = 
27rT2  = 


i(M) 


-z(l,i)log2:|f 

-i{l,i)z\\  (47) 


The  following  results  are  the  integrals  of  the  special  Green’s  function  terms,  using  the  same 
notation  as  above  for  the  terms  in  c;  the  conjugate  functions  which  must  be  combined  with  the 
terms  below  are  of  the  same  form  with  c  substituted  for  c.  The  coefficient  (1,2)  is  also  conjugated 
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=  [(^  “  c)  J(2,  c)  -  z  +  y/ z^  -  +  I{c)  log(2  +  -  a^)] 


(49) 


9c 


-z(l,z) 


Zy/^ 


y/c^  —  —  a?  cz  —  o?  I  \/c^  — 


4-2 


+  i(l,-i)[---(2,c)] 


47rri 

47rr2 


-i(l,i) 


J(z,c)  +c 


aj(^,c)  a/(c) 

3c  3c 


log(z  +  yj -  a^) 


+  i(l,-z)[---(z,c)] 


(50) 


8.4  Test  Problem  for  Antiplane  Shear  Loaded  Crack 

The  test  problem  is  the  antiplane  shear  fracture  mechanics  problem.  The  antiplane  shear  problem  is 
governed  by  Laplace’s  equation  V^it;(x,  y)  for  the  equilibrium  equation  in  terms  of  the  displacement 
function  w{x,y).  The  two  shear  stresses  (denoted  by  Txz  and  Ty^  where  z  is  normal  to  the  plane) 
correspond  to  the  gradients  of  the  potential  in  the  x,  y  directions.  Thus,  the  formulated  Green’s 
function  corresponds  to  the  exact  solution  for  a  large  plate  containing  a  central  crack  remotely 
loaded  by  the  applied  shear  traction  d<pjdn.  Equilibrium  of  the  boundary  tractions  is  given  by  the 
side-condition 

£^...0  (5.) 

which  must  be  satisfied  for  the  solution  to  the  BIE  to  exist. 

The  test  problem  is  a  square  plate  containing  a  central  crack  of  length  2a.  The  shear  modulus 
is  taken  to  be  unity.  The  selected  BEM  model  for  this  problem  is  given  in  Figure  1.  The  mesh 
consists  of  twenty-four  boundary  elements,  as  shown.  The  crack  is  taken  to  have  an  angle  (3  relative 
to  the  global  x-axis.  The  test  problem  is  loaded  by  shear  values  of  4-5  units  on  the  top  surface,  -5 
units  on  the  bottom  surface,  and  zero  on  the  side  surfaces. 

The  validation  solutions  will  be  the  shear  stress  intensity  factor,  referred  to  as  the  Mode  III  stress 
intensity  factor  in  fracture  analysis.  The  intensity  factor  for  the  infinite  plate  loaded  by  a  constant 
Tyz  =  To  at  infinity  is  given  by  the  following  result  [9]. 

Km  =  ToV^ra  (52) 

In  the  case  of  the  angled  crack  the  solution  is  found  in  [10]  to  be  given  as  follows. 

Km  =  To  cos  I3y/im  (53) 
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The  exact  displacement  field  solution  for  points  on  the  crack  surfaces  (corresponding  to  the  values 
of  the  potential  in  the  field)  is  given  in  [9]  as 


^  ^  I  _  /x\2 
a  fjLy  \aJ 

The  singular  stress  terms  near  the  crack  tip  are  given  as  follows 
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(54) 
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The  BEM  validation  solutions  are  for  a  normalized  total  plate  width  of  W/2a  =  100.  The  solution 
for  the  stress  intensity  factor  is  judged  to  be  equivalent  to  that  of  the  infinite  crack.  The  BEM  code 
reproduces  the  exact  solution  for  the  stress  intensity  factor,  the  gradient  along  the  crack,  and  the 
displacement  of  the  crack.  The  accuracy  relative  to  the  above  analytical  solutions  for  the  infinite 
body  is  0.004%,  and  can  be  further  decreased  by  increasing  W/2a.  The  BEM  code  is  therefore  taken 
to  be  validated. 


9  Application  to  the  Electromagnetic  Problem 

The  target  problem  is  the  prediction  of  the  gradient  of  the  three  dimensional  magnetic  field  created 
by  steady  current  flow  through  a  two  dimensional  plate  containing  one  or  two  cracks  coming  from  a 
hole.  The  electrical  field  problem  in  two  dimensions  is  solved  using  the  Green’s  function  formulation 
of  the  previous  sections.  We  then  apply  three  dimensional  potential  theory  to  get  the  induced 
magnetic  field  normal  to  the  plate.  All  three  dimensional  aspects  of  the  current  flow  problem  in  the 
finite  thickness  plate  are  ignored  and  the  two  dimensional  current  flow  is  integrated  over  the  finite 
thickness  of  the  plate  to  get  the  three  dimensional  magnetic  vector  at  a  specified  field  location.  The 
normal  component  of  the  magnetic  field  is  detected  at  that  location  using  a  fine-scale  magnetometer 
type  of  device. 

The  following  few  equations  summarize  the  application  of  the  special  Green’s  function  to  the 
problem  of  the  electromagnetic  (EM)  field  problem  of  steady  current  flow  through  a  body  with  a 
crack.  The  surfaces  will  be  assumed  to  be  perfectly  insulating  except  for  the  surfaces  of  current 
injection  and  removal.  The  conservation  of  charge  is  a  side  condition  for  this  Neumann  problem. 
The  following  derivations  are  based  on  the  equations  and  definitions  in  [8]. 

The  vector  electrical  field  for  the  steady-electrical  conduction  problem  E{q)  is  taken  to  be  the 
negative  gradient  of  a  harmonic  function. 


E{q)  =  -VV{q)  (56) 

The  current  vector  J{q)  is  scaled  from  the  electrical  field  by  the  conductivity  cr  which  is  taken  to  be 
a  constant  for  the  body. 

J{q)  =  aE(q)  (57) 

The  steady  current  flow  problem  means  that  the  electrostatic  potential  V{q)  is  harmonic,  —  0. 

The  boundary  conditions  for  the  conductor  are  given  in  terms  of  the  current  flow  density  across 
the  boundary  defined  by  its  outward  normal  n  unit  vector. 

-  dV 

J{Q)-n{Q)  =  -a—  =  crg{Q)  (58) 
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For  insulated  portions  of  the  boundary,  g{Q)  =  0.  Thus,  we  have  a  Neumann  problem  to  solve  with 
the  BIE,  for  the  potential  V,  The  BIE  for  this  problem  is  given  below, 

j\y{Q)  -  =  J^^{P,Q-,a)g{Q)dS{Q)  (59) 

The  physical  interpretation  of  the  requirement  for  the  solvability  of  this  Fredholm  equation  of  the 
First  Kind  is  that  the  total  current  flow  is  conserved.  When  we  account  for  the  fact  that  parts  of 
the  regular  boundary  are  also  insulated,  we  get  the  flnal  form  of  the  EM-BIE  for  our  application. 

ljV{Q)-V{P)]^^^^dSiQ)  =  '^{P,Q-,a)g{Q)dSiQ)  (60) 

The  solution  of  the  EM-BIE  is  by  the  use  of  boundary  element  modeling  of  the  surface  and  the 
imposition  of  boundary  data  interpolatives  for  each  boundary  element.  The  EM-BIE  code  is  based 
on  the  algorithm  in  [6]  except  that  we  now  use  quadratic  variations  for  all  boundary  data,  straight- 
line  boundary  elements,  and  exact  integration  of  the  boundary  element  terms. 

The  vector  magnetic  fleld  B  at  an  arbitrary  point  in  space  but  caused  by  the  steady  current  is 
derived  according  to  the  law  of  Biot-Savart  [8],  We  will  take  the  three  dimensional  field  point  for 
the  determination  of  the  magnetic  field  as  the  point  c(x).  The  magnetic  field  is  given  as  the  curl 
of  the  current  field  divided  by  the  distance  between  the  coil  point  and  the  integration  point.  The 
result,  written  for  the  three  dimensional  problem,  is  given  as  follows. 


where  /zq  is  the  free  space  permeability  constant.  Note  that  the  gradient  is  taken  at  the  free-space 
point  c  and  not  the  integration  point  q.  The  gradient  could  be  taken  outside  of  the  integral  and  lead 
to  a  different  potential  formulation  for  the  non-steady  state  solution.  We  will  leave  the  gradient 
inside  the  integration  for  our  developments. 

The  result  can  also  be  written  in  terms  of  the  electrical  current  potential  V{q),  which  substitution 
is  valid  for  the  steady-state  electromagnetic  formulation,  as  given  in  Eq.  56. 


At  this  time  we  can  expand  the  gradient  operator  as  follows 


Vc  X 


\  r(c,  q)  J  ^  y  r{c,  q)  J  r(c,  q) 


(63) 


However,  the  last  term  above  is  zero  due  to  taking  the  curl  of  a  gradient  of  a  scalar  field.  We  can 
then  write  the  following  form  in  terms  of  the  derivatives  at  q. 


Vc  X 


\  »’(c,9)  / 


(64) 


The  derivatives  are  now  all  at  the  point  q  and  we  can  drop  the  notation  of  the  derivative  point. 
We  now  use  Stokes’  theorem  to  convert  the  volume  integral  into  the  equivalent  surface  integral,  as 
follows.  The  form  of  Stoke’s  theorem  for  the  current  application  is  now  given. 


j  V  xFdV  =  jnxFdS 


(65) 
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where  u  is  the  surface  normal.  The  model  problem  is  a  plate  with  contours  defined  in  the  x~y  plane 
and  two  parallel  surfaces  in  the  2:-direction.  Thus,  if  we  take  the  2:-component  of  the  left-hand  side 
of  Eq.  65,  the  resulting  surface  integrals  are  on  the  two-dimensional  contours  for  the  plate  being 
analyzed.  It  is  this  component  of  the  gravity  field  that  we  wish  to  model.  The  second  identity 
converts  the  cross-product  term  to  the  following  dot-product  term. 


63  *  (n  X  F)  =  ((§3  X  n)  •  F 

=  i’F 


(66) 


and  thus 


^  63  •  (V  X  F)dV  =  ^  i-FdS 


(67) 


Substituting  the  results  from  Eq.  62  and  63  into  Eq.  67,  we  obtain  the  necessary  boundary 
integrals  for  the  B3  component  of  the  magnetic  field  in  terms  of  the  current  flow  along  the  boundaries 
of  the  two-dimensional  plate. 


The  modeling  approach  that  is  to  be  used  for  the  EM-applications  is  to  be  quasi-three  dimen¬ 
sional.  We  take  the  current  to  be  defined  for  two  dimensional  slices  parallel  to  the  surface  of  the 
plate.  Each  slice  has  a  two  dimensional  surface  and  crack  geometry,  although  the  crack  geometry 
can  be  changed  for  each  layer  to  simulate  a  three  dimensional  crack.  The  quasi-three  dimensional 
approach  neglects  flow  components  in  the  direction  normal  to  the  surface  of  the  plate.  This  approx¬ 
imation  may  be  valid  if  changes  in  the  geometry  in  the  direction  normal  to  the  plane  of  the  plate 
are  small. 

The  tangent  vector  in  the  above  magnetic  field  problem  is  now  unique  and  is  the  vector  tangent  to 
the  two  dimensional  surface  which  has  been  employed  in  the  EM-BIE  formulation.  The  tangential 
derivatives  of  the  potential  V  at  the  boundary  can  therefore  be  determined  from  the  EM-BIE 
solution,  including  points  on  the  crack  tip.  Analytical  derivation  of  the  special  Green’s  function  to 
get  the  tangential  derivatives  is  used  to  compute  the  singular  current  flow  at  the  crack  tip.  However, 
the  above  form  is  still  in  three  dimensional  form  and  its  integration  will  be  discussed  in  the  next 
section. 


%Vi,) .  -fJjViQ)  -  ±  I 


dVjQ)  dnP,Q-,a) 


dn 


dt 


dSiQ)  (69) 


The  magnetic  field  at  c(x)  is  derived  from  a  volume  potential  converted  to  the  surface  tangential 
current  flow  problem.  If  the  current  flow  was  truly  two  dimensional,  the  two  dimensional  result  for 
the  magnetic  field  would  convert  the  above  equation  to  the  log(2  -  c)  form  used  for  the  fundamental 
solution.  Clearly,  the  quasi-three  dimensional  approach  contains  the  most  critical  three  dimensional 
elements  of  the  magnetic  field  problem  and  can  probably  serve  well  in  the  current  application,  is 
spite  of  the  current  flow  approximations  used. 


9.1  The  Biot-Savart  Integral 

The  Biot-Savart  law  [11]  is  used  to  calculate  the  z-component  of  the  magnetic  field  from  the  tangential 
derivative  of  the  scalar  potential,  V.  The  integral  from  Eq.  68  is  written  again  in  the  local  coordinate 
system  of  the  measurement  device,  as  follows 


B,{x)  =  - 


/xotr/i 
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W(Q)-t 
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dS{Q) 


(70) 
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pickup  coil 


Figure  2:  Integration  geometry  between  pickup  coil  and  boundary  segment 


where  h  is  the  thickness  of  the  plate.  The  BEM  algorithm  calculates  dV/dt  as  a  piecewise  linear 
result  since  we  are  taking  V(Q)  to  be  represented  as  piecewise  quadratic.  Letting  Vt  =  VV  'i  be  the 
tangential  component  of  the  gradient  of  the  scalar  potential,  we  have  for  the  line  segment  along  the 
arc  length  s,  where  the  surface  is  taken  to  have  unit  thickness  {dS  =  1  •  ds). 

Vt{Q)  =  bi-^b2S  (71) 

This  integration  is  done  for  all  boundary  elements  except  for  the  crack  surface  F.  Substituting 
this  Vt  into  the  integral  for  Bz{c)  gives  the  contribution  to  the  magnetic  field  at  point  c  due  to  all 
boundaries  except  for  the  crack  surface,  F. 
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To  carry  out  the  integration,  the  geometrical  relations  for  r  and  s  must  be  determined. 
From  Figure  2  we  can  see  that 

D 

r  =  - 

cos  6 

D  =  y/z^+fP' 

,  rd9  DdO 

(XS  -  rr:  - 

COS  6  cos'^Q 

s  =  Dtan^l^^  =  D(tan0  -  tan^i) 

Therefore,  for  the  boundaries  other  than  the  crack,  the  integral  is  now 
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(74) 


where  Bi  =  bi  —  b2Dtaxi0i. 

Since  the  surface  5  is  divided  into  n  segments,  the  summation  of  this  equation  ~  evaluated  over 
each  segment  for  which  Vt  is  defined  by  a  single  linear  equation  (each  defined  by  a  different  61  and 
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62)  -  will  determine  Bz  at  the  pickup  coil 
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The  angles  61  and  02  correspond  to  the  angles  between  the  perpendicular  D  and  each  endpoint  of 
the  segment  as  the  summation  moves  in  a  positive  sense  around  the  boundary  (material  on  left). 

For  points  on  the  crack  surface,  the  singular  behavior  of  Vt  at  the  crack  tip  requires  special 
consideration.  The  complete  analytic  representation  of  Vt  for  crack  problems  is  given  by  [12] 


/i(^) 

\/x^  — 


+  f2{x) 


(76) 


where  /i(x)  and  f2{x)  are  analytic  functions  of  the  real  position  x  on  the  crack  and  a  is  the  crack 
half-length. 

The  first  term  in  Eq.  76  contains  the  discontinuity  at  the  crack  tip  while  the  second  term  is 
needed  to  match  far  field  (away  from  the  crack)  boundary  conditions.  On  the  crack  surface,  Eq.  76 
can  be  rewritten  as 


Vilr  = 


+  f2{x) 


y/a  —  Xy/a  +  X 
■^iFi{x)  +  f2ix) 


(77) 


where  R  is  the  distance  measured  from  the  crack  tip  at  x  =  a.  For  the  region  near  the  positive  crack 
tip,  the  behavior  of  Vt  is  known  to  follow 


Vi  = 


K+  .  (3 

\ - sin  -7 

2 


(78) 


where  K'^  is  the  potential  intensity  factor  (PIF)  of  the  crack  on  the  positive  x-axis  side.  Note  that 
at  =  0  (crack  tip),  Vt  =  oo. 

If  there  were  a  crack  on  the  opposite  side  (negative  x-direction)  of  the  rivet,  then  there  would 
be  a  corresponding  Vt  term  proportional  to  K~ ,  The  sin  function  multiplier  reflects  a  geometrical 
dependence  of  Vt  on  /?,  the  angle  measured  from  the  crack  tip  axis  to  the  evaluation  point  (see 
Figure  1).  For  points  on  the  upper  surface  of  the  crack,  /3  =  tt  and  the  multiplier  is  a  (+1).  For 
points  on  the  lower  surface  of  the  crack,  /?  =  — tt  and  the  multiplier  is  a  (-1).  Therefore,  combining 
Eq.  77  and  Eq.  78,  the  value  (real  part  only)  of  the  term  multiplying  the  singular  term,  Fi(x),  can 
be  expressed  in  terms  of  the  PIF 

Fi(a)  =  ±^-/2(a)  (79) 

v27r 

Due  to  the  singular  behavior  of  Vt  at  the  crack  tip,  numerical  integration  will  be  necessary  to 
evaluate  the  magnetic  field  contribution  of  the  crack.  Values  of  Vt  will  be  evaluated  for  the  upper  and 
lower  crack  surfaces  using  the  BIE  program  and  then  the  magnetic  field  will  be  calculated  through 
a  numerical  integration  of  the  previously  stated  Biot-Savart  relation  (Eq.  70). 


9.2  Numerical  integration  of  the  Biot-Savart  integral 

Numerical  integration  over  the  crack  is  accomplished  using  discretization  of  the  crack  boundary  into 
elements.  By  utilizing  a  coordinate  transformation,  Vt  can  be  expressed  in  terms  of  nodal  values 
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regular  element 

Figure  3:  Crack  divided  into  regular  and  quarter  point  (crack-tip)  elements 


Figure  4:  Quadratic  interpolation  functions  for  mapping  the  regular  and  crack  tip  elements  into  the 
intrinsic  ^-space 


and  interpolation  functions  (or  shape  functions)  [13]  of  an  intrinsic  coordinate  Once  mapped 
into  ^-space,  Gaussian  quadrature  numerical  integration  can  be  used  to  evaluate  the  integral.  The 
present  formulation  will  use  quadratic  interpolation  functions  Nm{0  which  requires  three  nodes  per 
boundary  element. 

The  crack  will  be  modeled  using  one  crack-tip  element  containing  the  crack-tip  node  and  one 
regular  element,  each  element  being  defined  by  three  nodes  (see  Figure  3).  The  regular  element 
is  a  straight-forward  use  of  the  interpolation  functions  but  the  crack-tip  element  will  require  some 
modification  to  accommodate  the  singularity  at  the  crack  tip.  The  crack- tip  element  contains  the 
crack  tip  node  as  its  first  node  and  the  mapping  uses  the  quarter-point  formulation  [13]  to  account 
for  the  inverse  square-root  singular  behavior  at  the  crack  tip.  The  other  regular  element  will  use  the 
general  form  of  the  shape  functions  to  map  them  to  ^-space. 

The  general  form  of  the  quadratic  interpolation  functions  used  to  map  the  function  to  be  inte¬ 
grated  into  the  ^-space  are  as  follows; 

iVi(0  = 

N2{0  =  1-^"  (80) 

NziO  =  ^^(^  +  1) 

The  parent  element  along  the  ^-axis  with  uniformly  spaced  nodes  at  ^  —  1,  0,+l  and  the  parent 

shape  functions  are  shown  in  Figure  4. 

The  spatial  coordinate  transformation,  s  =  s(^)  maps  the  parent  element  ranging  from  f  —  1 
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to  ^  =  4-1  onto  each  of  the  crack  boundary  elements.  A  set  of  three  parent  shape  functions  are 
created  for  each  element.  It  is  desirable  that  this  mapping  be  one-to-one  so  that  each  point  in  the 
parent  element  maps  onto  one  point  in  the  real  element.  In  particular,  the  boundary  nodes  of  the 
parent  element,  ^  =  1,  must  map  onto  the  boundary  nodes  of  the  real  element  to  ensure  continuity 
at  the  interelement  boundaries.  For  the  quadratic  element,  the  mapping  takes  the  form: 

3 

S  =  (81) 

fc=l 

where  si^S2,  and  S3  are  the  corresponding  real-space  nodal  values  of  the  variable  being  mapped  from 
the  real  element  to  the  parent  element.  It  can  be  seen  that  this  equation  maps  the  parent  nodes 
onto  the  real  nodes: 


^  =  -1  — >  s  =  Si 

^  0^s  =  S2  (82) 

^  =  -hi  ^  s  =  S3 


Only  the  values  of  si,S2,  and  S3  for  each  element  need  to  be  changed  for  each  mapping.  The 
coordinate  transformation  into  ^-space  of  the  integral  results  in  the  following  sum 

[  F{s)ds  ->  T  r  N^{OFkJ{s.Od^  (83) 

where  Fk  is  the  value  of  the  function  at  the  node  and  J(s,0  the  Jacobian  of  the  mapping  of 
the  real  surface  onto  the  reference  surface. 


J{s.O 


ds 

di 


dNk 


(84) 


The  Jacobian  can  be  written  out  explicitly  for  use  with  quadratic  shape  function  as  follows. 

Jis,  0  =  (e  -  0  -  2^52  +  (^  +  isa)  (85) 

For  the  regular  element,  the  location  of  the  second  node  is  midway  between  the  endpoint  nodes  of 
the  element,  S2  =  {si  4-  s^)/2,  so  the  resulting  mapping  is 

S  =  -\m-Osi  +  {l-e)s2  +  \ai  +  Os3  (86) 

_  (1  ~  0^1  +  (1  +  0^3 

2 


with  the  corresponding  Jacobian 


j  _  ds  _  S3  —  Si 
~d^~  2 


(87) 


For  Si  =  0  and  S3  =  L,  the  contribution  to  Bz  of  the  regular  element  on  the  upper  crack  surface 
becomes 


(88) 
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+  A^3(0 


S2 


+N2{i) 


j{s,m 

di 


This  result  can  now  be  numerically  integrated  using  regular  Gaussian  quadrature  techniques  [13]. 

For  the  crack-tip  element  on  the  upper  crack  surface,  the  standard  mapping  to  {-space  needs 
to  be  modified  to  accommodate  the  l/\/s  behavior  of  Vi(s)  at  the  crack  tip.  In  what  follows,  the 
crack  tip  point  will  be  mapped  at  s  =  0.  The  general  approach  will  be  to  split  the  Vt{s)  term  into 
a  singular  term  multiplied  times  a  non-singular  coefficient.  The  form  of  Vt(5)  follows  from  Eq.  77. 


Vt{s)\r 


y/s^/2L  —  s 


(89) 


To  map  Vt(5)  into  {-space,  the  singular  term  will  use  an  inverse  mapping  relation  and  the  non¬ 
singular  term  will  use  the  quadratic  shape  functions.  To  determine  the  inverse  mapping  relation, 
we  start  with  the  l/\/s  behavior  which  can  be  represented  by  placing  the  midnode  of  the  quadratic 
element  at  the  quarter  point  location 

S2  =  \{SZ  -  5a)  (90) 

Therefore,  choosing  si  =  0,  S2  =  S3  =  L  in  Eq.  86,  the  quarter  point  mapping  is  given  by 
the  following  result 

s  =  (l-e')|  +  ^^(l  +  0i  (91) 

which  has  a  corresponding  inverse  mapping  given  as  follows. 


c  = 


(92) 


The  Jacobian  J  is  given  by 


-2^|  +  i(l  +  20I 

f(l  +  0 


(93) 


Note  that  at  the  s  =  0({  —  —1)  point  the  Jacobian  also  equals  zero.  This  characteristic  is 
important  when  evaluating  the  integral  of  the  mapped  function.  The  l/y/s  term  is  represented  in 
{-space  by  solving  for  l/\/s  in  terms  of  {  in  Eq.  92 


1-  ^  2 
v/5  ”  ^/I(?  + 1) 


(94) 


Therfore,  the  Bz  contribution  due  to  the  crack-tip  element  on  the  upper  crack  surface  is  as  follows 


(95) 
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Note  that  as  ^  ^  -1  the  singular  term  (^  +  1  in  denominator)  goes  to  infinity  but  is  canceled 
by  the  Jacobian  term  (^  +  1  in  numerator)  going  to  zero,  thereby  making  the  overall  function  well 
behaved  in  the  mapped  space.  All  values  of  F,  r,  and  s  are  known  except  for  F  at  the  crack  tip  node 
which  has  the  singularity.  But  at  the  crack  tip,  the  value  of  the  function  is  the  PIF(A'+) 


FM  =  ^  (96) 

All  boundaries  are  now  taken  care  of  and  the  following  section  summarizes  the  experimental  valida¬ 
tion  of  the  measurement  model  based  on  this  development. 


9.3  Experimental  Validation 

A  measurement  model  has  been  developed  using  this  BIE/EM  formulation.  The  program  simulates 
the  scanning  of  a  Superconducting  Quantum  Interference  Device  (SQUID)  magnetometer  [14]  over  a 
sample  containing  a  crack.  The  SQUID  magnetometer  uses  small  pickup  coils  to  sense  magnetic  field 
above  the  sample.  Injecting  a  uniform  DC-current  in  the  sample  causes  the  current  to  be  parallel  to 
the  sample’s  surface  under  the  pickup  coil.  The  associated  magnetic  field  is  mostly  tangential  to  the 
plate  surface  for  scans  located  centrally  and  with  small  lift-off  distances.  A  flaw  in  the  specimen  will 
perturb  the  current  distribution  and  produce  a  vertical  magnetic  field  component  that  can  then  be 
detected  by  the  SQUID.  When  the  SQUID  is  scanned  two-dimensionally  over  the  sample,  a  magnetic 
field  map  is  produced,  revealing  a  flaw  signature  that  commonly  has  a  dipolar  shape  (maximum  and 
minimum  peaks). 

As  a  measurement  reference,  electro  discharge  machined  (EDM)  slots  and  saw  cuts  are  used  in 
calibrating  NDE  systems  for  detection  of  cracks.  A  combination  of  a  drilled  hole  with  an  EDM  slot 
is  an  approximation  to  a  crack  emanating  from  an  aircraft  fastener  hole.  Fabrication  of  test  samples 
made  this  way  is  simple  and  controllable  making  it  easy  to  build  a  test  set  representing  the  range 
of  conditions  that  are  of  interest.  But  measurements  with  NDE  instruments  [15]  have  shown  that 
the  instrument  response  is  not  necessarily  the  same  as  that  from  a  fatigue  crack  of  the  same  size 
and  geometry.  It  is  possible  that  crack  closure  may  cause  electrical  conductivity  across  parts  of  the 
fatigue  crack  thus  causing  a  reduced  signal  response  over  an  electrically  insulated  slot.  The  BIE 
formulation  used  in  the  SQUID  measurement  model  is  for  a  closed  crack  (no  width)  but  electrically 
insulated  along  the  entire  length.  These  issues  must  be  taken  into  consideration  when  validating 
the  model  against  experimental  measurements  on  open  cracks  (i.e.,  slots)  and  actual  fatigue  cracks 
which  may  have  closure.  For  this  work,  magnetic  field  map  shape  comparisons  are  the  beisis  for 
validation  and  provides  most  of  the  information,  except  amplitude,  needed  to  accomplish  this. 

Following  are  descriptions  of  two  experimental  measurements  on  fabricated  samples  representa¬ 
tive  of  those  commonly  used  in  SQUID  NDE. 

Case  I:  The  first  validation  experiment  used  a  5  mA  DC-current  injected  into  a  75  mm  x  150  mm 
X  0.03  mm  copper  clad  circuit  board  containing  a  15  mm  x  0.03  mm  slot  cut  with  a  scalpel  (see 
Figure  5).  Although  this  setup  does  not  provide  completely  uniform  current  injection  across  the 
sample  (transverse  to  the  slot)  due  to  the  point  source  electrodes,  the  region  around  the  centrally 
located  slot  should  be  relatively  uniform. 


no 


75X 150X  0.03  mm  Cu  plate 


Figure  5:  Case  I:  Experimental  setup  for  DC-current  injection  in  copper  plate  with  slot 


Figure  6  shows  the  contour  map  resulting  from  a  2-D  scan  over  the  sample  with  the  lines  AA\ 
BB\  and  CC’,  Profiles  along  these  lines  on  this  map  are  compared  to  those  calculated  by  the 
measurement  model.  The  scaling  factor  between  the  measured  signal  response  (a  voltage)  and  the 
calculated  magnetic  flux  was  determined  using  these  profiles.  Magnetic  field  shape  characteristics 
(matching  of  peaks  and  valleys)  will  reflect  the  accuracy  of  the  measurement  model  since  a  difference 
between  the  model  and  experiment  would  show  up  as  a  mismatch  of  the  profiles  at  either  the  edge 
or  the  crack  locations. 

The  measurement  model  is  in  very  good  agreement  with  experiment  for  all  profiles.  A  slight 
mismatch  at  the  right  side  of  AA’  can  be  seen  and  is  most  likely  due  to  a  small  variation  in  liftoff 
during  a  scan  (i.e.  sample  not  level). 

Case  II:  The  second  validation  measurement  was  similar  to  the  first  with  a  5  mA  DC-current 
injected  into  a  100  mm  x  150  mm  x  0.03  mm  copper  clad  circuit  board  containing  a  9  mm  diameter 
hole  with  a  9  mm  x  0.03  mm  slot  on  one  side  of  the  hole  (see  Figure  7).  However,  the  SQUID 
used  for  this  measurement  was  a  high-Tc  system  [16]  which  used  a  planar  gradiometer  (pickup  and 
balance  coils  in  same  plane)  instead  of  axially  along  a:.  The  scans  were  taken  in  an  unshielded 
environment. 

Figure  8  shows  the  contour  map  resulting  from  a  2-D  scan  over  the  sample.  The  quadrupolar 
shape  results  from  the  planar  gradiometer  acting  as  a  spatial  differentiater  and,  by  taking  the 
derivative  of  the  dipole-shaped  magnetic  signal,  results  in  a  quadrupolar  shape.  Since  the  planar 
gradiometer  was  oriented  parallel  to  the  plate  edges,  the  edge  signal  is  approximately  zero  and 
therefore,  only  that  part  of  the  scan  above  the  hole/slot  was  used  in  these  comparisons.  Profiles 
along  lines  on  this  map  are  compared  to  those  calculated  by  the  measurement  model. 

Figure  8  shows  the  profile  comparisons,  after  scaling,  corresponding  to  the  lines  AA\  BB’,  and 
CC’  on  the  contour  map.  The  measurement  model  is  in  very  good  agreement  with  experiment  for 
all  profiles.  Again,  a  difference  between  the  model  and  experiment  would  show  up  as  a  mismatch  of 
the  profiles  at  the  location  of  the  hole/crack. 


10  Summary  and  Conclusions 

The  problem  formulation  begins  with  a  boundary-integral  representation  for  Laplace’s  equation  in 
the  plane.  The  presence  of  a  single,  straight  crack  in  the  plane  is  represented  by  a  special  Green’s 
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Figure  6:  Case  I:  Magnetic  field  contour  map  and  profiles  comparing  experiment  with  the  BEM 
model 


100X150X0.03  mmCu  plate 


Figure  7:  Case  II:  Experimental  setup  for  DC-current  injection  in  copper  plate  with  hole/slot 
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Figure  8:  Case  II:  Magnetic  field  contour  map  and  profiles  comparing  experiment  with  the  BEM 
model 


function  corresponding  to  zero  flux  across  the  crack.  The  Green’s  function  used  provides  a  very 
effective  BIE  model  for  the  zero-flux  boundary  condition  as  the  crack  then  is  a  part  of  the  integral 
equation  kernels  and  is  not  a  part  of  the  modeled  boundary.  A  boundary  element  implementation 
for  the  harmonic  BIE  uses  quadratic  data  variation  on  piecewise  linear  boundary  elements.  The 
BEM  model  with  the  Green’s  function  is  validated  for  the  antiplane  shear  crack  problem  for  which 
exact  solutions  exist. 

The  steady-state  electromagnetic  problem  suitable  for  calculating  the  magnetic  field  for  a  cracked 
plate  with  current  injection  has  been  formulated  in  terms  of  the  derived  BIE.  The  gradient  of  the 
potential  yields  the  electrical  field  vector.  The  electromagnetic  field  equations  yield  the  magnetic 
field  at  arbitrary  three-dimensional  space  points,  using  the  law  of  Biot-Savart.  The  law  of  Biot-Savart 
result  is  given  in  terms  of  two  boundary  integrals  of  current  flow  tangent  to  the  boundaries  which 
must  be  evaluated  for  all  boundary  points,  including  the  crack.  The  cracked  plate  is  represented 
as  a  finite  thickness  plate  with  a  two-dimensional  field  so  that  the  BEM  results  can  be  used  in  the 
three-dimensional  magnetic  field  computation. 

The  magnetic  field  integrals  on  the  boundary  points  are  formulated  as  Gaussian  quadratures 
for  the  quadratic  BEM  data.  The  crack  surface  requires  a  special  integration  algorithm  due  to 
the  singular  conditions  at  each  crack  tip.  Isoparametric  boundary  element  algorithms  are  used 
for  the  magnetic  field  boundary  integrals.  The  crack  tip  element  uses  the  quarter-point  mapping 
developed  for  BEM  traction  models  of  fracture  problems  which  contain  the  leading  term  in  the 
singular  behavior.  The  mapping  is  shown  to  contain  the  leading  term  exactly. 

The  application  problem  is  the  magnetic  field  for  the  steady-state  current  flow  through  two- 
dimensional  plate  with  a  crack.  The  magnetic  field  is  the  normal  component  of  the  field,  which  is 
the  component  detected  by  an  experimental  SQUID  magnetometer.  Experimental  data  was  obtained 
for  two  cracked  plates.  The  first  is  a  copper  plate  with  an  EDM  slot  representing  a  crack.  The  second 
is  a  copper  plate  with  a  circular  hole  and  a  single  edge  crack  normal  to  the  hole  surface. 

The  new  BEM  results  for  the  magnetic  field  are  used  to  simulate  a  SQUID  scan  over  the  plate, 
along  three  scan  lines  parallel  and  normal  to  the  current  flow  field.  The  agreement  between  the 
BEM  results  and  the  experimental  data  is  excellent. 
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APPENDIX  B: 


C-PROGRAM 


/*  PROGRAM:  biescan.c  */ 

/*  Tony  Ewing,  Nov  96  */ 

/*  Program  to  solve  the  magnetic  forward  problem  using  BIE,  */ 

/*  generating  the  simulated  magnetic  field  data  due  to  circular  rivet 
hole  */ 

with  crack  which  would  be  detected  by  a  single-coil  magnetometer.  */ 

#include  <stdlib.h> 

#include  <stdio.h> 

#include  <math.h> 

#define  pi  3.141592653589793 
#define  Nh  200 
#define  Nc  200 

main ( ) 

{ 

/*  scan  */ 

char  scan_in[20] ,hole_in[20] , crack_in[20]  ; 
int  i, j , k, 1, widsteps, Ingsteps, njmp, cnum; 

double  scnwid, trcewid, scnlng, scninc, xstart, ystart, offset, xinc, theta; 
double  coilrad, coilrd2, area, loff, product , baseline, base; 
double  Inginc, widinc, X, Y, Z, calarea, Rcen; 
double  mu, cnst, rtd, Bmax, Bmin, Bpp; 

double  intflx [200]  [200] ; 

/*  hole  */ 

int  nc, nr, sn, eseg, hseg; 

double  dml, dm4 , kay, kayprm, dell, del2, snphil, csphil , phil , bl, blp, b2 ; 
double  dm21, dm2 2, sy, sx, snphi2, csphi2, phi2, s, seal, dm3 1, dm3 2, Bh, totf Ixhl; 
double  Bhtot , shatnx, shatny, rhol, rholx, rholy, rldtsnx,  rldtsny,  Dl,  rldtsn; 
double  rho2, rho2x, rho2y, r2dtsnx, r2dtsny, philp, phi2p, betal, beta2; 
double 

shattx, shatty, rldtstx, rldtsty, rldtst, r2dtstx, r2dtsty, r2dtst, signl, sign2 

/ 

double  a,b, c, d, e, fl, targl, targ2, xprint , yprint, xbegin, ybegin,phi, snphi; 
double  rl, r2, Dlprm, h, sig, calib, totflxh2 ; 

double  seg[Nh+l] [3] ; 
double  Vth [Nh+1] [3] ; 
double  flxh[200] [200] ; 
double  flxt [200] [200] ; 

/*  crack  */ 

int  gs, snj , sgint, cseg, cnod; 

double  nl, n2, n3, rxil, rxi2, rxi3, Vtxi, jl, j2, j3, rhonl, rhon2, rhon3, jxi, AC; 
double  sndl, snd2, snd3, rhonlx, rhonly, rhon2x, rhon2y, rhon3x, rhon3y; 

double  Vtxil, Vtxi2, Vtxi3, jxil, jxi2, jxi3, Lx, Ly, L, VtRl, VtR2, VtR3; 
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double  Be, Betot , totf Ixcl, totf lxc2 , holrad, biesign; 

double  node [Nc+1] [3] ; 
double  Vtc[Nc+l]; 
double  flxc[200]  [200] ; 
double  I5[6] ,xi5[6] ,w5[6] ; 

/*  print  */ 
int  ip, jp, print; 

/*  math  functions  */ 

double  cos (double) , sin (double) , acos (double) , tan (double) 
double  f abs (double) , pow (double, double) , sqrt (double) ; 
double  log (double) ; 

FILE  *fpin; 

FILE  ^fpinl; 

FILE  *fpin2; 

FILE  *fpout ; 

FILE  *fpoutl; 

FILE  *fpout2; 


/*  open  scanning  input  file  */ 

printf ( "Input  scan  input  file  name:  "); 
scanf (”%s", scan_in) ; 

if ( (fpin=fopen (scan_in, "r") ) ==NULL) 

{ 

printf ("cannot  open  file\n”); 
exit  (1) ; 

} 


/*  open  BIE  segments  input  file  for  hole  */ 

f scanf (fpin, ”%s”, hole_in) ; 

if  (  (fpinl=fopen  (hole_in,  "r")  )  —NULL) 

{ 

printf ( "cannot  open  file\n"); 
exit  ( 1 ) ; 

} 


/*  open  BIE  segments  input  file  for  crack  */ 

f scanf  (fpin,  ”%s",  crack___in)  ; 

if  (  (fpin2==fopen  (crack^in,  ”r")  )  ==NULL) 

{ 

printf ("cannot  open  file\n"); 
exit ( 1 ) ; 

} 

/*  open  output  file  */ 

if  (  (fpout==f open  ("scan. out",  "a")  )==NULL) 

{ 
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printf {"\nCannot  open  file\n\n") ; 
exit (1) ; 

} 

if  (  ( f pout l==f open  ("scanl. out",  "w")  )“NULL) 

{ 

printf ( "\nCannot  open  file\n\n"); 
exit (1) ; 

} 

if ( ( f pout 2=f open { "scan2 . out", "w") )==NULL) 

{ 

printf ("\nCannot  open  file\n\n"); 
exit (1) ; 

} 


/*  read  in  boundary  info,  SQUID  parameters  and  x-y  scan  geometry  */ 


fscanf (fpin, ”%lf", &holrad) ; 
fscanf {fpin, "%lf ", &AC) ; 
fscanf  (fpin,  "%lf  ",'&sig)  ; 
fscanf (fpin, "%lf ", &h) ; 
fscanf (fpin, "%lf", &calib) ; 
fscanf (fpin, "%d", &eseg) ; 
fscanf (fpin, "%d", &hseg) ; 
fscanf (fpin, "%d", &njmp) ; 
fscanf { fpin, "%d" , &cnod) ; 
fscanf (fpin, "%d", &cseg) ; 
fscanf (fpin, "%d", &sgint) ; 
fscanf (fpin,  "%lf ", Scoilrad) ; 
fscanf ( fpin,  "%lf ", &lof f ) ; 
fscanf (fpin,  "%lf ", &baseline) ; 
fscanf (fpin,  "%lf ", &scnwid) ; 
fscanf (fpin,  "%lf ", &scnlng) ; 
fscanf (fpin,  "%lf ", &lnginc) ; 
fscanf (fpin,  "%lf ", &widinc) ; 
fscanf (fpin,  "%lf ", &offset ) ; 


/*  read  in  x, y~coordinates  of  boundary  segment  endpoints  for  hole  */ 
/*  njmp  is  the  number  of  boundary  jumps  -  e,g.  2  closed  boundaries  is 
*/ 

/*  njmp=l  */ 

for (nr=l;nr<= (eseg+hseg+njmp+l) ;nr++) 

{ 

for  (nc=l;nc<=2;nc+-f) 

{ 

fscanf (fpinl, "%lf ", &seg [nr] [nc] ) ; 

/*  seg[nr] [nc] = (seg [nr] [nc] ) /lOO. ;*/ 

} 

for (nc=l;nc<=2;nc++) 


fscanf (fpinl, "%lf", &Vth[nr] [nc] ) ; 

printf ("Vth[%d] [%d]=  %f \n" , nr , nc, Vth [nr] [nc] ) ; 

} 

} 

/*  read  in  x, y-coordinates  of  nodes  for  crack  (3  per  crack  element)  */ 
/*  nseg  is  the  total  number  of  segments  (2  nodes  per  segment)  */ 
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/*  sgint  is  the  number  of  segments  per  element  -  eg  3  nodes/2 
segments/  */ 

/*  1  element.  */ 

print f ("cnod=  %d\n”, cnod) ; 

for {nr=l;nr<=cnod;nr++) 

{ 

for (nc=l;nc<=2;nc++) 

{ 

fscanf (fpin2, "%lf ”, &node [nr] [nc] ) ; 

} 

fscanf (fpin2, ”%lf ”, &Vtc [nr] ) ; 
print f ( ”Vtc [%d] =  %f \n”, nr, Vtc [nr] ) ; 

} 


/*  Define  constants  */ 


mu=pi*4 . Oe-07; 
rtd=180. /pi; 
cnst=mu/ (4 . 0*pi) ; 
coil rd2=coilrad* coil rad; 
calarea=pi*coilrd2; 

widsteps=scnwid/widinc;  /*  number  of  scan  traces  -  y  direction  */ 
lngsteps=scnlng/lnginc;  /*  number  of  steps  along  x-scan  */ 

/*  write  statements  for  ASCII  output  */ 

/*fprintf {fpout,  "%i\n", Ingsteps) ; 
fprintf { fpout,  "%i\n”, widsteps) ; */ 
fprintf ( fpout 1,  ”%i\n", Ingsteps) ; 
fprintf (fpoutl,  ”%i\n" , widsteps ) ; 
fprintf (fpout 2,  "%i\n”, Ingsteps) ; 
fprintf (fpout 2,  "%i\n", widsteps) ; 

/^***************^^******  START  SCAN  ***★*********★★*****★*/ 

for (i=0;i<=widsteps;i++)  /*  ’loop  to  step  along  y  */ 

{ 

printf("loop  %d  of  %d\n”, i, widsteps) ; 
xstart=-scnlng/2 . ; 

ystart=-of f set- (scnwid/2 . ) +i* (widinc) ; 

for ( j=0; j<=lngsteps; j++)  /*  loop  to  step  along  x  */ 

xinc=xstart+ j  *lnginc; 
totf lxhl=0 . 0; 
totf Ixcl^O . 0; 
totflxh2=0 , 0; 
totflxc2=0 . 0; 
base=baseline; 

for (cnum=l; cnum<=2; cnum++)  /*  loop  for  gradiometer  */ 

{ 

if {cnum==2)  base=0 . 0; 

for (k=0; k<=6; k++)  /*  loop  for  coil  integration  -  Roth  */ 

{ 

theta= (pi/3 . ) *k; 
if (k==0) 

{ 

Rcen=0 . 0; 

area=0 . 25*calarea; 
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} 

else 

{ 

Rcen=0 . 81650*coilrad; 
area=0 . 125*calarea; 

} 

X=Rcen*cos (theta) +xinc; 

Y=Rcen*sin (theta) tystart; 

/***********★***********  hole  begin  ***********************^*/ 

/*  Calculate  vectors  between  source  and  field  points  */ 

Bhtot=0.0; 

for (sn=l; sn< (eseg+hseg+njmp+1 ) ; sn++) 

{ 

if (sn==eseg+l)  sn++; 

sy=seg[sn+l] [2]“seg[sn] [2] ; 

sx=seg[sn+l]  [l]’-seg[sn]  [1]; 

s=sqrt (pow(sx, 2) +pow{sy, 2) ) ; 

shatnx==sy/s; 

shatny=“sx/s; 

shattx=sx/s; 

shatty=sy/s; 

/*  If  modeling  rivot  hole  with  crack,  since  node  position  from  BIE  is 
in  */ 

/*  crack  coordinates,  need  to  shift  back  to  global  coordinates  by  */ 
AC  (crack  length)  -  only  on  X  now  -  subtract  holrad  from  rholx  and 
rho2x*/ 


if(cnod==0)  rholx=X-seg [sn] [1] ; 

else  rholx=X-seg[sn] [1] -holrad;  /^Ist  error:  subtract  hole 
radius  */ 

rholy=Y-seg [sn]  [2]; 

rhol=sqrt (pow(rholx, 2) +pow(rholy, 2) ) ; 

rl=sqrt (pow (rhol, 2) +pow ( (lof f+base)  ,  2 ) ) ; 

rldtsnx=rholx*shatnx; 

rldtsny=rholy*shatny; 

rldtstx=rholx*shattx; 

rldtsty=rholy*shatty; 

rldtst=rldtstx+rldtsty; 

rldtsn=rldtsnx+rldtsny; 

/*  Determine  sign  on  angle  using  dot  product  with  s-tangent  */ 
if  (rldtst==0 .0)  signl==-l .  0; 
else  signl=- (rldtst/fabs (rldtst) )  ; 

Dl=fabs (rldtsn) ; 

Dlprm=sqrt (pow (Dl, 2) +pow ( (loff+base) , 2) ) ; 
if (cnod==0)  rho2x=X-seg [sn+1] [1]; 

else  rho2x=X-seg [sn+1] [1] -holrad;  part  2  of  1st  error  */ 

rho2y=Y-seg [sn+1] [2] ; 

rho2=sqrt (pow (rho2x, 2) +pow (rho2y, 2) ) ; 

r2=sqrt (pow (rho2, 2) +pow ( (loff+base) ,2)); 

r2dtsnx=rho2x*shatnx; 

r2dtsny=rho2y*shatny; 

r2dtstx=rho2x*shattx; 

r2dtsty=rho2y*shatty; 
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r2dtst=r2dtstx+r2dtsty; 

if (r2dtst==0 .0)  sign2=-l . 0/ 
else  sign2=- (r2dtst/fabs (r2dtst) ) ; 

phil=signl*acos (Dlprm/rl) ; 
phi2=sign2*acos  (Dlprin/r2)  ; 


/*  Determine  linear  constants  from  tangential  derivatives  */ 

b2=((Vth[sn] [2] )-(Vth[sn] [l]))/s; 
bl= (Vth [sn] [1] ) ~b2*Dlprm*tan (phil) ; 
targl=tan ( (pi/4 . ) + {phil/2 . ) ) ; 
targ2=tan ( (pi/4 . ) + (phi2/2 . ) ) ; 

Bh= (bl*log (targ2/targl) )  +  (b2* {r2-rl)  )  ; 

Bhtot=Bh+Bhtot ; 

}  /*  end  loop  for  hole  bie  segments  */ 

if (cnum==l ) 

{ 

totf lxhl=Bhtot*area+totf Ixhl; 

} 

else 

{ 

totf lxh2=Bhtot*area+totflxh2 ; 

} 

/*★****★***********★★★***  hole  end  -crack  begin 

-A- -A- + T*r  / 

/*  Set  up  Gauss  points  and  corresponding  weights  arrays  for 
integration  */ 

xi5[0]=-sqrt { ( 35 . +2 . *sqrt ( 70 . ))/63.); 
xi5[l]=-sqrt( (35.-2. *sqrt (70. ) )/63.  )  ; 
xi5 [ 2 ] =0 . 0 ; 
xi5 [3] =-xi5 [1] ; 
xi5 [4]=-xi5 [0]  ; 

w5[0]=5103./{50.* (322.+13.*sqrt (70.  )  )  )  ; 
w5 [1]=5103./ (50.* (322.-13. *sqrt (70. ) ) ) ; 
w5[2]-128./225. ; 
w5[3]=w5[l] ; 
w5[4]=w5[0]  ; 


/*  Calculate  vectors  between  source  and  field  points  */ 
Bctot=0 . 0; 
snj=0; 

for (sn=0;sn<(cseg/2) ;sn++)  /*  2nd  error:  divide  by  2  */ 

{ 

if (sn%sgint==0)  snj=snj+l; 
if(sn==0)  snj=0; 

/*  5-pt  gaussian  quadrature  integration  loop  */ 
for {gs=0;gs<=4 ;gs++) 
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{ 

nl=.5*xi5 [gs] * (xi5 [gs]-l . )  ; 
n2=(l.-xi5 [gs] ) * (1 .+xi5 [gs] ) ; 
n3=. 5*xi5 [gs] * (xi5 [gs] +1. ) ; 

/*  Evaluate  terms  in  integral  */ 

rhonlx=X- (holrad+node [2*sn+l+snj  ]  [ 1] )  ; 
rhon2x=X- (holrad+node [2*sn+2+snj ] [1] ) ; 
rhon3x=X- (holrad+node [2*sn+3+snj ] [1] ) ; 
rhonly=Y-node[2*sn+l+snj] [2]; 
rhon2y=Y-node [2*sn+2+snj ]  [2]  ; 
rhon3y=Y-node [2*sn+3+snj ]  [2] ; 
rhonl=sqrt (pow(rhonlx, 2 . ) +pow (rhonly, 2 . ) ) ; 
rhon2=sqrt (pow (rhon2x, 2 . ) +pow(rhon2y, 2 . ) ) ; 
rhon3=sqrt (pow(rhon3x, 2 . ) +pow{rhon3y, 2.)); 

rxil= (1. / (sqrt (pow ( (loff+base) ,2.) +pow(rhonl, 2 .)))); 
rxi2= (1 . / (sqrt (pow ( (loff+base) , 2. ) +pow(rhon2, 2. ) ) ) ) ; 
rxi3= (1. / (sqrt (pow( (loff+base) ,2. ) +pow (rhon3, 2 . ) ) ) ) ; 

Vtxil=nl*Vtc[2*sn+l+snj] ; 

Vtxi2=n2*Vtc[2*sn+2+snj]; 

Vtxi3=n3*Vtc [2*sn+3+snj ] ; 

/*  Determine  sign  change  due  to  numerical  integration  direction  */ 

/*  starting  at  craclc  tip  nodes  rather  than  bie  integration  being  */ 
/*  positive  direction  with  material  always  on  left  */ 
Lx=node[2*sn+3+snj] [1] -node [2*sn+l+snj ] [1] ; 

Ly=node [2*sn+3+snj ] [2] -node [2*sn+l+snj ] [2] ; 

L=sqrt (pow (Lx, 2) +pow (Ly, 2 ) ) ; 

if (node [2*sn+l+snj ] [2] >0.0)  biesign=1.0; 
else  biesign=-l . 0; 

/*  Check  to  see  if  on  a  crack-tip  element  or  regular  element  */ 
if (sn%sgint==0) 

{ 

VtR2=sqrt (AC-fabs (node [2*sn+2+snj ] [1])); 
VtR3=sqrt(AC-fabs(node[2*sn+3+snj] [1] ) )  ; 
jxi=biesign*sqrt (L) ; 

15 [gs] = (Vtxil*rxil+Vtxi2*rxi2*VtR2+Vtxi3*rxi3*VtR3) * jxi; 

} 

else 

{ 

15 [gs] = (Vtxil*rxil+Vtxi2*rxi2+Vtxi3*rxi3) *biesign* {L/2 . ) ; 

} 


}  /*  end  gaussian  quadrature  loop  */ 

Bc=w5[0]*I5[0]+w5[l]*I5[l]+w5[2]*I5[2]+w5[3]*I5[3]+w5[4]*I5[4] ; 
/*  if (ystart==0)  printf ("Crack  loop  [%d] ,  Bc=  %e\n", sn, Be) ; */ 

Bctot=Bc+Bctot ; 

}  /*  end  loop  for  crack  bie  segments  */ 
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if (cnum==l) 

{ 

totf lxcl=Bctot*area+totflxcl; 

} 

else 

{ 

totflxc2=Bctot*area+totflxc2; 

} 

}  /*  end  loop  for  coil  integration  */ 

}  /*  end  loop  for  gradiometer  */ 

flxh[i] [j]=cnst*h*sig*calib* (totflxh2~totflxhl) ; 
flxc[i] [j]=cnsf^h*sig*calib* ( tot f lxc2-totflxcl ) ; 
flxt [i] [j]=flxh[i] [j]+flxc[i] [j] ; 

}  /*  end  loop  for  scan  along  x  */ 

)  /*  end  loop  for  increment  along  y  */ 

/***★*★*★******★*★★★****★**  *****^*****************^^^^*^^ 
/*  Set  initial  values  within  output  range  */ 

Bmax=flxt [1] [1] ; 

Bmin=flxt [1]  [1]  ; 


/*  Loop  to  print  output  file  */ 
for (ip=0;ip<= (widsteps) ;ip++) 

{ 

xbegin=-scnlng/2 . ; 

yprint=-of f set- {scnwid/2 . ) tip* (widinc) ; 
for ( jp=0; jp<=lngsteps; jp++) 

{ 

/*  Find  max/min  of  scan  */ 

if (flxt [ip] [ jp] >Bmax)  Bmax=flxt [ip] [ jp] ; 

if (flxt [ip] [jp]<Bmin)  Bmin=f Ixt [ip] [ jp] ; 


xprint=xbegin+jp*lnginc; 

/*  fprintf (fpout,  "%f  %f  %e\n",  xprint , yprint , flxh [ip] [ jp] ) ; */ 

/*  fprintf (fpoutl,  "%f  %f  %e\n”, 

xprint, yprint, intflx [ip] [ jp] ) ; */ 

fprintf (fpout2,  "%f  %f  %e\n”,  xprint , yprint , f Ixt [ ip] [ jp] ) ; 

/*  fprintf (fpout2,  "%e\n", flxt [ip] [ jp] ) ; */ 

} 

} 

Bpp=Bmax-Bmin ; 

fprintf  (fpout,  **%e\n",  Bpp)  ; 


f close (fpin) ; 
fclose (fpinl) ; 
fclose (fpin2) ; 
fclose (fpout ) ; 
fclose (fpoutl) ; 
fclose ( fpout2 ) ; 
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}  /*  end  of  roainO  */ 
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APPENDIX  C: 


FLUXA^OLTAGE  CALIBRATION  OF  SQUID  MAGNETOMETERS 


FluxA^oltage  Calibration  of  SQUID  Magnetometers 

The  SQUID  magnetometer  outputs  a  voltage,  F,  proportional  to  the  magnetic  flux, 
0,  through  the  pickup  coil.  Many  SQUID  systems  use  gradiometers  consisting  of  a 
pickup  coil  connected  serially  to  a  counter-wound  balance  coil  (Figure  C.l).  This  is  done 
to  reject  noise  from  distant  field  sources  (uniform)  since  the  flux  through  each  coil  would 
be  approximately  the  same  resulting  in  zero  net  flux.  For  closer  sources,  such  as 
scanning  a  sample,  the  flux  in  stronger  in  the  pickup  coil  than  in  the  compensation  coil 
resulting  in  a  net  flux  and  hence,  an  output  voltage. 


balance  coil 


pick-up  coil 


Figure  C.  1  Gradiometer  with  pick-up  coil  and  counter-wound  balance  coil 


To  use  SQUIDs  as  accurate  quantitative  measuring  tools,  correct  instrument 
calibration  is  necessary  to  determine  the  absolute  magnetic  field  being  measured.  One 
approach  is  to  scan  the  SQUID  over  a  standard  source  for  which  the  magnetic  field 
spatial  distribution  can  be  determined  analytically  so  that  comparisons  can  be  made.  A 
commonly  used  source  is  a  current  carrying  wire  (see  Figure  C.l).  From  Ampere’s  Law, 
it  can  be  determined  that  the  vertical  component  of  magnetic  field,  B2,  in  the  x-direction  a 
distance  z  above  a  current  (I)  carrying  wire  is  given  by 
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(C.l) 


Figure  C.l  Calibration  using  a  current  carrying  wire 

Differentiating  this  relation  gives  the  location  of  the  maximum  and  minimum  peaks  of  Bz 
resulting  in  relationship  for  the  peak-to-peak  separation  (see  Figure  C.3) 

dB,  z^-x^ 

dx  (x^  +  z^J  (C.2) 

bx-lz. 


AX 

Figure  C.3  Dipolar  signal  showing  peak-to-peak  separation  and  amplitude 


So,  the  peak-to-peak  field  value,  Bpp^  is  determined  by  the  difference  of  at  x-  ±z 
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This  can  be  used  to  estimate  the  vertical  location  (z)  of  the  pickup  coil  above  the  wire 
through  a  measurement  of  Bpp.  This  is  only  an  approximation  since  it  assumes  a  single 
point  detector  above  an  infinite  wire.  To  model  the  output  from  a  gradiometer  above  a 
finite  wire,  the  magnetic  field  must  be  integrated  over  the  areas  of  the  pickup  and  balance 
coils. 

Integration  schemes  can  be  found  in  [21]  and  used  to  integrate  the  magnetic  field 
from  a  finite  length  wire,  length  21  carrying  current  I 


It  is  usually  difficult  to  know  z  accurately  since  the  pickup  coils  are  located  inside 
the  dewar.  To  address  the  uncertainties  in  the  system  variables,  a  procedure.  Calibration 
Requiring  Approximate  Parameters  [22],  has  been  developed  to  determine  the  factor 
scaling  the  output  voltage  to  magnetic  flux.  The  procedure  uses  optimization  [23]  over 
liftoff,  coil  tilt,  gradiometer  baseline,  coil  radius,  and  the  calibration  factor  to  “match”  the 
experimentally  measured  voltage  map  to  the  magnetic  flux  map  calculated  using  Eq. 
(C.4).  The  resulting  voltage-to-flux  calibration  factor  is  independent  of  magnetic  field 
source  and  liftoff  distance  and  therefore,  can  be  used  for  all  measurements  using  that 
system. 
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So  far  in  this  discussion,  the  calibration  factor  has  scaled  voltage  to  net  magnetic 
flux  {0-  units  of  weber).  It  is  common  practice  to  state  this  calibration  factor  as  a 
voltage  to  magnetic yzcW  {B  —  units  of  tesla)  number.  This  is  done  by  converting  the  net 
flux  back  to  an  approximate  magnetic  field  value  by  dividing  by  the  area  of  the  pickup 
coil.  This  hov^ever  introduces  error  into  the  resulting  value  since  this  conversion  must 
assume  a  uniform  magnetic  field  across  the  coil  area  and  ignores  the  balance  coil 
contribution.  This  is  useful  for  those  interested  in  approximate  magnetic  field  values  that 
are  not  dependent  on  the  system’s  coil  geometry,  as  flux  is. 


single  coil  -  uniform  B  gradiometer  -  nonuniform  B 


<D=  ^B-dA 

AO=  {B^-dA^ 

J 

-  P! 

coil 

pickup  coil 

hafance  coil 

=>B  =  — 

=>B  =  1 

A 

These  issues  affect  only  quantitative  amplitude  values  since  the  result  is  an 
optimized  calibration  factor  used  to  scale  the  theoretical  result  to  the  experimental 
measurement.  The  shape  of  the  magnetic  maps  is  unchanged  by  this  process  and  does  not 
affect  qualitative  analysis  of  the  images. 
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